From 593581aee2fbb40ad0f95718c67590b0885dd901 Mon Sep 17 00:00:00 2001 From: Elod Csirmaz Date: Tue, 5 Mar 2024 00:11:55 +0000 Subject: [PATCH] Add PathTube and tube-like surfaces --- README.md | 10 +++ openscad_py.py | 189 ++++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 188 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 6153210..4b1b546 100644 --- a/README.md +++ b/README.md @@ -23,3 +23,13 @@ becomes translate(v=[0, 0, 1]) { cube(size=[1, 1, 2], center=false); } ``` +## Notable convenience functions + +- Usual computational geometry functions on the `Point` class that work in an arbitrary number of dimensions. Overloads for algebraic operators + +- `Cylinder.from_ends` constructs a cylinder between two given points in space + +- `Polyhedron.tube` creates a tube-like or toroid polyhedron from a 2D array of points + +- `PathTube` creates a tube-like or toroid polyhedron from an arbitrary path + diff --git a/openscad_py.py b/openscad_py.py index 9b2cdb0..64c3f86 100644 --- a/openscad_py.py +++ b/openscad_py.py @@ -1,5 +1,6 @@ from typing import Union as TUnion +from typing import List import math import numpy as np @@ -89,12 +90,16 @@ class Point: def allclose(self, p: 'Point') -> bool: return self.c.shape == p.c.shape and np.allclose(self.c, p.c) - def angle(self, p: 'Point') -> float: - """Return the angle between two vectors, in degrees""" + def angle(self, p: 'Point', mode: str = "deg") -> float: + """Return the angle between two vectors in degrees or radians""" r = self.dot(p) r = r / self.length() / p.length() r = math.acos(r) - return r / math.pi * 180. + if mode == "rad": + return r + if mode == "deg": + return r / math.pi * 180. + raise ValueError("Unknown mode") def rotate(self, coords, angle: float) -> 'Point': """Rotate. coords is a list of 2 coordinate indices that we rotate""" @@ -181,6 +186,7 @@ class Header: self.quality = quality def render(self): + # See https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Other_Language_Features#Circle_resolution:_$fa,_$fs,_and_$fn if self.quality == 'draft': return "" if self.quality == 'mid': @@ -191,7 +197,9 @@ class Header: class Cube(Object): - """A 3D primitive, cube""" + """A 3D primitive, cube. + Creates a cube in the first octant. When center is true, the cube is centered on the origin.""" + # https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/The_OpenSCAD_Language#cube def __init__(self, size: TUnion[list, Point], center: bool = False): self.size = Point.c(size) @@ -202,7 +210,9 @@ class Cube(Object): class Sphere(Object): - """A 3D primitive, sphere""" + """A 3D primitive, sphere. + Creates a sphere at the origin of the coordinate system.""" + # https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/The_OpenSCAD_Language#sphere def __init__(self, r): self.r = r @@ -213,7 +223,9 @@ class Sphere(Object): class Cylinder(Object): - """A 3D primitive, cylinder""" + """A 3D primitive, cylinder. + Creates a cylinder or cone centered about the z axis. When center is true, it is also centered vertically along the z axis.""" + # https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/The_OpenSCAD_Language#cylinder def __init__(self, h, r=None, r1=None, r2=None, center: bool = False): self.height = h @@ -234,7 +246,7 @@ class Cylinder(Object): length = v.length() assert length != 0 z = Point([0, 0, 1]) - r = z.cross(v).norm() + r = z.cross(v) rangle = v.angle(z) if r.length() == 0: # The cylinder is in the Z direction @@ -242,13 +254,168 @@ class Cylinder(Object): p1 = p2 rangle = 0 r = z + else: + r = r.norm() return cls(h=length, r=radius, center=False).rotate(a=rangle, v=r).move(p1) -class Circle(Object): - """A 2D primitive, circle""" +class Polyhedron(Object): + """A 3D primitive, a polyhedron defined by a list of points and faces.""" + # See https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/The_OpenSCAD_Language#polyhedron + # Nonplanar faces should be triangulated by opensCAD + + def __init__(self, points: List[TUnion[list, Point]], faces: List[list], convexity: int = 10): + self.points = [Point.c(p) for p in points] + self.faces = faces + self.convexity = convexity - def __init__(self, r, fn=None): + @classmethod + def tube(cls, points: List[List[TUnion[list, Point]]], convexity: int = 10): + """Construct a tube-like polyhedron from a 2D array of points. + Each row of points must be oriented clockwise when looking at the pipe at the start inwards. + The rows of points form loops. + """ + rows = len(points) + row_len = len(points[0]) + point_list = [] + point_map = {} # { (row_ix,col_ix) -> list_ix, ... + for row_ix, row in enumerate(points): + for col_ix, point in enumerate(row): + point_map[(row_ix, col_ix)] = len(point_list) + point_list.append(point) + + faces = [] + + # Side faces + for row_ix in range(1, rows): + for col_ix in range(1, row_len): + faces.append([ + point_map[(row_ix, col_ix-1)], + point_map[(row_ix, col_ix)], + point_map[(row_ix-1, col_ix)], + point_map[(row_ix-1, col_ix-1)] + ]) + faces.append([ + point_map[(row_ix, row_len-1)], + point_map[(row_ix, 0)], + point_map[(row_ix-1, 0)], + point_map[(row_ix-1, row_len-1)] + ]) + + # Starting cap + faces.append([point_map[(0,x)] for x in range(row_len)]) + + # Ending cap + faces.append([point_map[(rows-1,row_len-1-x)] for x in range(row_len)]) + + return cls(points=point_list, faces=faces, convexity=convexity) + + def render(self) -> str: + faces_list = [f"[{','.join([str(x) for x in face])}]" for face in self.faces] + return f"polyhedron(points=[{','.join([p.render() for p in self.points])}], faces=[{','.join(faces_list)}], convexity={self.convexity});" + + +class PathTube(Object): + """Creates a tube-like or toroid polyhedron from a path (list of points).""" + + def __init__(self, points: List[TUnion[list, Point]], radius: float, fn: int, convexity: int = 10): + self.points = [Point.c(p) for p in points] + self.radius = radius + self.fn = fn # number of sides + self.convexity = convexity + + def process(self, debug: bool = False) -> Polyhedron: + points_rows = [] + + for ix, point in enumerate(self.points): + if debug: print(f"//LOOP {ix}: {point.render()}") + + if ix == 0: + # Start of the path + v = self.points[1].sub(point) # vector toward the first point + z_point = Point([0,0,1]) + seam = v.cross(z_point) # Track a seam along the pipe using this vector pointing from the middle line + if seam.length() == 0: # v is in the z direction + seam = Point([1,0,0]) + seam = seam.norm() + seam2 = v.cross(seam).norm() + if debug: print(f"//Start. v={v.render()} seam={seam.render()} seam2={seam2.render()}") + points = [] + for i in range(self.fn): + a = math.pi*2*i/self.fn + points.append((seam*math.cos(a) + seam2*math.sin(a))*self.radius + point) + points_rows.append(points) + if debug: print(f"// Row: {', '.join([p.render() for p in points])}") + + elif ix == len(self.points) - 1: + # End of the path + v = point.sub(self.points[-2]) + seam2 = v.cross(seam).norm() + if debug: print(f"//End. v={v.render()} seam={seam.render()} seam2={seam2.render()}") + points = [] + for i in range(self.fn): + a = math.pi*2*i/self.fn + points.append((seam*math.cos(a) + seam2*math.sin(a))*self.radius + point) + points_rows.append(points) + if debug: print(f"// Row: {', '.join([p.render() for p in points])}") + + else: + # Middle of the path + # (p[-1]) -va-> (p[0]) -vb-> (p[1]) + va = point.sub(self.points[ix-1]).norm() # vector incoming to this elbow + vb = self.points[ix+1].sub(point).norm() # vector going out from this elbow + if debug: print(f"//Middle. va={va.render()} vb={vb.render()}") + # Get the vector perpendicular to va that points to the inside of the cylinder around va according + # to the elbow at p[0]. This is the component of vb in a basis defined by va. + vdot = va.dot(vb) + vb_proj = va.scale(vdot) # The projection of vb onto va + vb_perp = vb.sub(vb_proj) # This is perpendicular to va + if debug: print(f"// vb_proj={vb_proj.render()} vb_perp={vb_perp.render()}") + va_inner = vb_perp.norm() + + va_proj = vb.scale(vdot) + va_perp = va.sub(va_proj) + if debug: print(f"// va_proj={va_proj.render()} va_perp={va_perp.render()}") + vb_inner = va_perp.scale(-1).norm() # Here we want to project -va onto vb + if debug: print(f"// va_inner={va_inner.render()} vb_inner={vb_inner.render()}") + + # The new seam on vb (seam_b) has the same angle to vb_inner as it had on va to va_inner + seam_angle = seam.angle(va_inner, mode="rad") + # need to figure out the sign of the angle + if seam_angle != 0: + if va_inner.cross(seam).dot(va) < 0: + seam_angle = -seam_angle + vb_inner2 = vb.cross(vb_inner).norm() + seam_b = vb_inner*math.cos(seam_angle) + vb_inner2*math.sin(seam_angle) + if debug: print(f"// seam={seam.render()} seam_b={seam_b.render()}") + + vangle = va.scale(-1).angle(vb, mode="rad") + long_inner = (vb-va).norm().scale(1/math.sin(vangle/2)) + # long_inner is the long axis of the elliptic intersection between the cylinders around va and vb + short = va.cross(long_inner).norm() # the short axis of the ellipse + if debug: print(f"// long_inner={long_inner.render()} short={short.render()} vangle={vangle/math.pi*180}(deg) seam_angle={seam_angle/math.pi*180}(deg)") + points = [] + for i in range(self.fn): + # We draw the ellipse according to long_inner and short, but use seam_angle to get the right points + a = math.pi*2*i/self.fn + seam_angle + points.append((long_inner*math.cos(a) + short*math.sin(a))*self.radius + point) + points_rows.append(points) + if debug: print(f"// Row: {', '.join([p.render() for p in points])}") + + seam = seam_b + + return Polyhedron.tube(points=points_rows, convexity=self.convexity) + + def render(self) -> str: + return self.process().render() + + +class Circle(Object): + """A 2D primitive, circle. + Creates a circle (or regular polygon) at the origin.""" + # https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/The_OpenSCAD_Language#circle + + def __init__(self, r: float, fn: TUnion[int, None] = None): self.r = r self.fn = fn # $fa, $fs, $fn @@ -279,7 +446,7 @@ class Polygon(Object): def render(self) -> str: return f"polygon(points=[{','.join([p.render() for p in self.points])}], convexity={self.convexity});" -# TODO polyhedron(points=[[],], faces[[p,],], convexity=) + # TODO https://docs.python.org/3/reference/datamodel.html#emulating-numeric-types