GNAL/scad/libraries/path_extrude.scad

196 lines
7.6 KiB
OpenSCAD
Raw Normal View History

// path_extrude.scad -- Extrude a path in 3D space
// usage: add "use <path_extrude.scad>;" to the top of your OpenSCAD source code
// Copyright (C) 2014-2019 David Eccles (gringer) <bioinformatics@gringene.org>
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
// Determine the projection of a point on a plane centered at c1 with normal n1
function project(p, c, n) =
p - (n * (p - c)) * n / (n * n);
// determine angle between two points with a given normal orientation
// see https://stackoverflow.com/questions/14066933/
// direct-way-of-computing-clockwise-angle-between-2-vectors
// dot = p1 * p2;
// det = (p1[0]*p2[1]*n1[2] + p2[0]*n1[1]*p1[2] + n1[0]*p1[1]*p2[2]) -
// (n1[0]*p2[1]*p1[2] + p1[0]*n1[1]*p2[2] + p2[0]*p1[1]*n1[2]);
// atan2(det, dot);
// determine angle between two planar points and a centre
// with a given normal orientation
function getPlanarAngle(p1, p2, c1, n1) =
let(p1 = p1-c1, n1=n1 / norm(n1), p2=p2-c1)
atan2((p1[0]*p2[1]*n1[2] + p2[0]*n1[1]*p1[2] + n1[0]*p1[1]*p2[2]) -
(n1[0]*p2[1]*p1[2] + p1[0]*n1[1]*p2[2] + p2[0]*p1[1]*n1[2]), p1 * p2);
function c3D(tPoints) =
(len(tPoints[0]) == undef) ? // single point
c3D([tPoints])[0] :
(len(tPoints[0]) < 3) ? // collection of 2D points
tPoints * [[1,0,0],[0,1,0]] :
tPoints; // 3D points
// translate a point (or points)
function myTranslate(ofs, points, acc = []) =
(len(points[0]) == undef) ?
myTranslate(ofs, [points])[0] :
[ for(i = [0:(len(points) - 1)])
[ for(d = [0:(len(points[0])-1)]) (ofs[d] + points[i][d])]];
// rotate a point (or points)
function myRotate(rotVec, points) =
let(rotX = [[1, 0, 0],
[0, cos(rotVec[0]), -sin(rotVec[0])],
[0, sin(rotVec[0]), cos(rotVec[0])]],
rotY = [[ cos(rotVec[1]), 0,-sin(rotVec[1])],
[ 0, 1, 0],
[ sin(rotVec[1]), 0, cos(rotVec[1])]],
rotZ = [[ cos(rotVec[2]), sin(rotVec[2]), 0],
[ sin(rotVec[2]), -cos(rotVec[2]), 0],
[0, 0, 1]])
(len(points[0]) == undef) ?
myRotate(rotVec, [points])[0] :
c3D(points) * rotX * rotY * rotZ;
// Determine spherical rotation for cartesian coordinates
function rToS(pt) =
[-acos((pt[2]) / norm(pt)),
0,
-atan2(pt[0],pt[1])];
function calcPreRot(p1, p2, p3) =
let(n1=p2-p1, // normal between the two points (i.e. the plane that the polygon sits on)
n2=p3-p2,
rt1=rToS(n1),
rt2=rToS(n2),
pj1=(p2 + myRotate(rt2, [[1e42,0,0]])[0]),
pj2=project(p=(p1 + myRotate(rt1, [[1e42,0,0]])[0]), c=p2, n=n2))
getPlanarAngle(p1=pj1, p2=pj2, c1=p2, n1=n2);
function cumSum(x, res=[]) =
(len(x) == len(res)) ? concat([0], res) :
(len(res) == 0) ? cumSum(x=x, res=[x[0]]) :
cumSum(x=x, res=concat(res, [x[len(res)] + res[len(res)-1]]));
// Create extrusion side panels for one polygon segment as triangles.
// Note: panels are not necessarily be planar due to path twists
function makeSides(shs, pts, ofs=0) =
concat(
[for(i=[0:(shs-1)]) [i+ofs, ((i+1) % shs + ofs + shs) % (shs * pts),
(i+1) % shs + ofs]],
[for(i=[0:(shs-1)]) [((i+1) % shs + ofs + shs) % (shs * pts),
i+ofs, (i + ofs + shs) % (shs * pts)]]);
// Concatenate the contents of the outermost list
function flatten(A, acc = [], aDone = 0) =
(aDone >= len(A)) ? acc :
flatten(A, acc=concat(acc, A[aDone]), aDone = aDone + 1);
// Linearly interpolate between two shapes
function makeTween(shape1, shape2, t) =
(t == 0) ? shape1 :
(t == 1) ? shape2 :
[for (i=[0:(len(shape1)-1)])
(shape1[i]*(1-t) + shape2[i % len(shape2)]*(t))];
// Extrude a 2D shape through a 3D path
// Note: merge has two effects:
// 1) Removes end caps
// 2) Adjusts the rotation of each path point
// so that the end and start match up
module path_extrude(exPath, exShape, exShape2=[],
exRots = [0], exScale = [1], merge=false, preRotate=true){
exShapeTween = (len(exShape2) == 0) ?
exShape : exShape2;
shs = len(exShape); // shs: shape size
pts = len(exPath); // pts: path size
exPathX = (merge) ? concat(exPath, [exPath[0], exPath[1]]) :
concat(exPath,
[exPath[pts-1] + (exPath[pts-1] - exPath[pts-2]),
exPath[pts-1] + 2*(exPath[pts-1] - exPath[pts-2])]);
exScaleX = (len(exScale) == len(exPath)) ? exScale :
[for (i = [0:(pts-1)]) exScale[i % len(exScale)]];
preRots = [for(i = [0:(pts-1)])
preRotate ?
calcPreRot(p1=exPathX[i], // "current" point on the path
p2=exPathX[(i+1)], // "next" point on the path
p3=exPathX[(i+2)]) :
0 ];
cumPreRots = cumSum(preRots);
seDiff = cumPreRots[len(cumPreRots)-1]; // rotation difference (start - end)
// rotation adjustment to get start to look like end
seAdj = -seDiff / (len(cumPreRots));
adjPreRots = (!merge) ? cumPreRots :
[for(i = [0:(pts-1)]) (cumPreRots[i] + seAdj * i)];
adjExRots = (len(exRots) == 1) ?
[for(i = [0:(len(adjPreRots)-1)]) (adjPreRots[i] + exRots[0])] :
[for(i = [0:(len(adjPreRots)-1)]) (adjPreRots[i] + exRots[i % len(exRots)])];
phPoints = flatten([
for(i = [0:(pts-1)])
let(p1=exPathX[i],
p2=exPathX[(i+1)],
n1=p2-p1, // normal between the two points
rt1=rToS(n1))
myTranslate(p1, myRotate(rt1, myRotate([0,0,-adjExRots[i]],
c3D(makeTween(exShape, exShapeTween, i / (pts-1)) *
exScaleX[i]))))
]);
if(merge){ // just the surface, no end caps
polyhedron(points=phPoints,
faces=flatten([
for(i = [0:(pts-1)])
makeSides(shs, pts, ofs=shs*i)
])
);
} else {
polyhedron(points=phPoints,
faces=concat(
flatten([
for(i = [0:(pts-2)])
makeSides(shs, pts, ofs=shs*i)
]),
concat( // add in start / end covers
[[for(i= [0:(shs-1)]) i]],
[[for(i= [(len(phPoints)-1):-1:(len(phPoints)-shs)]) i]]
)
));
}
}
myPathTrefoil = [ for(t = [0:(360 / 101):359]) [ // trefoil knot
5*(.41*cos(t) - .18*sin(t) - .83*cos(2*t) - .83*sin(2*t) -
.11*cos(3*t) + .27*sin(3*t)),
5*(.36*cos(t) + .27*sin(t) - 1.13*cos(2*t) + .30*sin(2*t) +
.11*cos(3*t) - .27*sin(3*t)),
5*(.45*sin(t) - .30*cos(2*t) +1.13*sin(2*t) -
.11*cos(3*t) + .27*sin(3*t))] ];
myPointsOctagon =
let(ofs1=15)
[ for(t = [0:(360/8):359])
((t==90)?1:2) * [cos(t+ofs1),sin(t+ofs1)]];
myPointsChunkOctagon =
let(ofs1=15)
[ for(t = [0:(360/8):359])
((t==90)?0.4:1.9) *
[cos((t * 135/360 + 45)+ofs1+45)+0.5,sin((t * 135/360 + 45)+ofs1+45)]];
//myPoints = [ for(t = [0:(360/8):359]) 2 * [cos(t+45),sin(t+45)]];
pts=[2,0,0.5];
/*translate([0,0,0]) {
path_extrude(exRots = [$t*360], exShape=myPointsOctagon,
exPath=myPathTrefoil, merge=true);
}*/