Add PathTube and tube-like surfaces

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