update RELEASE

README.md

@@ -1,8 +1,8 @@
-# dotSCAD 1.1
+# dotSCAD 1.3
 
-> Helpful modules and functions when playing OpenSCAD.
+> Reduce the burden of 3D modeling in mathematics. Compatible with OpenSCAD 2015.03 or laters.
 
-![dotSCAD](GoldenTaiwan.JPG)
+![dotSCAD](WhirlingTaiwan.JPG)
 
 [![license/LGPL](https://img.shields.io/badge/license-LGPL-blue.svg)](https://github.com/JustinSDK/lib-openscad/blob/master/LICENSE)
 
@@ -27,6 +27,8 @@ If OpenSCAD generates "WARNING: Ignoring unknown xxx function" or "WARNING: Igno
 
 Too many dependencies? Because OpenSCAD doesn't provide namespace management, I personally think that exposing dependencies is better than hiding them. In this way, users can have their own way to manage dependencies. How to categorize dependencies is up to you. For example, you can include your commonly-used modules and functions in "commonly_used.scad" and then `include <commonly_used.scad>;` in the ".scad" file of your project.
 
+If you really don't want to care about dependencies, `include <dotSCAD.scad>;` or `use <dotSCAD.scad>;` come to save you. 
+
 ## Documentation
 
 - 2D
@@ -39,6 +41,7 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [hexagons](https://openhome.cc/eGossip/OpenSCAD/lib-hexagons.html)
 	- [polytransversals](https://openhome.cc/eGossip/OpenSCAD/lib-polytransversals.html)
     - [multi_line_text](https://openhome.cc/eGossip/OpenSCAD/lib-multi_line_text.html)
+	- [voronoi2d](https://openhome.cc/eGossip/OpenSCAD/lib-voronoi2d.html)
 
 - 3D
 	- [rounded_cube](https://openhome.cc/eGossip/OpenSCAD/lib-rounded_cube.html)
@@ -49,6 +52,8 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [hull_polyline3d](https://openhome.cc/eGossip/OpenSCAD/lib-hull_polyline3d.html)
 	- [function_grapher](https://openhome.cc/eGossip/OpenSCAD/lib-function_grapher.html)
 	- [polysections](https://openhome.cc/eGossip/OpenSCAD/lib-polysections.html)
+	- [starburst](https://openhome.cc/eGossip/OpenSCAD/lib-starburst.html)
+	- [voronoi3d](https://openhome.cc/eGossip/OpenSCAD/lib-voronoi3d.html)
 	
 - Transformation
     - [along_with](https://openhome.cc/eGossip/OpenSCAD/lib-along_with.html)
@@ -63,6 +68,13 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [parse_number](https://openhome.cc/eGossip/OpenSCAD/lib-parse_number.html)
 	- [cross_sections](https://openhome.cc/eGossip/OpenSCAD/lib-cross_sections.html)
 	- [paths2sections](https://openhome.cc/eGossip/OpenSCAD/lib-paths2sections.html)
+	- [path_scaling_sections](https://openhome.cc/eGossip/OpenSCAD/lib-path_scaling_sections.html)
+	- [bijection_offset](https://openhome.cc/eGossip/OpenSCAD/lib-bijection_offset.html)
+	- [in_polyline](https://openhome.cc/eGossip/OpenSCAD/lib-in_polyline.html)
+	- [in_shape](https://openhome.cc/eGossip/OpenSCAD/lib-in_shape.html)
+	- [midpt_smooth](https://openhome.cc/eGossip/OpenSCAD/lib-midpt_smooth.html)
+	- [trim_shape](https://openhome.cc/eGossip/OpenSCAD/lib-trim_shape.html)
+	- [triangulate](https://openhome.cc/eGossip/OpenSCAD/lib-triangulate.html)
 	
 - Path
     - [arc_path](https://openhome.cc/eGossip/OpenSCAD/lib-arc_path.html)
@@ -74,6 +86,7 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
     - [golden_spiral](https://openhome.cc/eGossip/OpenSCAD/lib-golden_spiral.html)
     - [archimedean_spiral](https://openhome.cc/eGossip/OpenSCAD/lib-archimedean_spiral.html)
     - [sphere_spiral](https://openhome.cc/eGossip/OpenSCAD/lib-sphere_spiral.html)
+	- [torus_knot](https://openhome.cc/eGossip/OpenSCAD/lib-torus_knot.html)
 
 - Extrusion
     - [box_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-box_extrude.html)
@@ -96,6 +109,7 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [shape_path_extend](https://openhome.cc/eGossip/OpenSCAD/lib-shape_path_extend.html)		
 
 - 2D Shape Extrusion
+	- [bend_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-bend_extrude.html)
 	- [path_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-path_extrude.html)
 	- [ring_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-ring_extrude.html)
 	- [helix_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-helix_extrude.html)

RELEASE.md

@@ -1,3 +1,47 @@
+> Version numbers are based on [Semantic Versioning](https://semver.org/).
+
+# v1.3.2
+- All-in-one source file.
+  - You can use `include <dotSCAD.scad>;` or `use <dotSCAD.scad>;` if you really don't want to care about dependencies.
+
+- Bugfixes
+  - `along_with`: Wrong variable scope.
+  
+# v1.3.1
+- Bugfixes
+  - `in_polyline`: Wrong parameter name.
+  - `in_shape`: Missing dependency.
+  - `along_with`: Avoid warning when using 2D points.
+  
+# v1.3
+- New modules:
+  - [bend_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-bend_extrude.html)
+  - [voronoi2d](https://openhome.cc/eGossip/OpenSCAD/lib-voronoi2d.html)
+  - [voronoi3d](https://openhome.cc/eGossip/OpenSCAD/lib-voronoi3d.html)
+
+- New functions:
+  - [in_shape](https://openhome.cc/eGossip/OpenSCAD/lib-in_shape.html)
+  - [in_polyline](https://openhome.cc/eGossip/OpenSCAD/lib-in_polyline.html)
+  - [midpt_smooth](https://openhome.cc/eGossip/OpenSCAD/lib-midpt_smooth.html)
+  - [trim_shape](https://openhome.cc/eGossip/OpenSCAD/lib-trim_shape.html)
+  - [triangulate](https://openhome.cc/eGossip/OpenSCAD/lib-triangulate.html)
+
+- New parameters:
+  - `distance` of [shape_taiwan](https://openhome.cc/eGossip/OpenSCAD/lib-shape_taiwan.html)
+  - `epsilon` of [bijection_offset](https://openhome.cc/eGossip/OpenSCAD/lib-bijection_offset.html)
+  - `method` of [path_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-path_extrude.html)
+  - `method` of [along_with](https://openhome.cc/eGossip/OpenSCAD/lib-along_with.html)
+   
+# v1.2
+- New modules and functions:
+  - [starburst](https://openhome.cc/eGossip/OpenSCAD/lib-starburst.html)
+  - [torus_knot](https://openhome.cc/eGossip/OpenSCAD/lib-torus_knot.html)
+  - [bijection_offset](https://openhome.cc/eGossip/OpenSCAD/lib-bijection_offset.html)
+  - [path_scaling_sections](https://openhome.cc/eGossip/OpenSCAD/lib-path_scaling_sections.html)
+
+- Others
+  - Avoid warnings when using newer versions of OpenSCAD after 2015.03.
+
 # v1.1.1
 - Bugfixes
   - `m_rotation` returns an identity matrix if `a` is 0.
@@ -23,7 +67,6 @@
 - New Parameters:
   - added `v` parameter to [rotate_p](https://openhome.cc/eGossip/OpenSCAD/lib-rotate_p.html)
 
-
 - Improved Performance:
     - [path_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-path_extrude.html)
     - [align_with](https://openhome.cc/eGossip/OpenSCAD/lib-along_with.html)  

all/dotSCAD.scad

@@ -0,0 +1,79 @@
+include <along_with.scad>;
+include <arc.scad>;
+include <arc_path.scad>;
+include <archimedean_spiral.scad>;
+include <archimedean_spiral_extrude.scad>;
+include <bend.scad>;
+include <bend_extrude.scad>;
+include <bezier_curve.scad>;
+include <bezier_smooth.scad>;
+include <bezier_surface.scad>;
+include <bijection_offset.scad>;
+include <box_extrude.scad>;
+include <circle_path.scad>;
+include <cross_sections.scad>;
+include <crystal_ball.scad>;
+include <ellipse_extrude.scad>;
+include <function_grapher.scad>;
+include <golden_spiral.scad>;
+include <golden_spiral_extrude.scad>;
+include <helix.scad>;
+include <helix_extrude.scad>;
+include <hexagons.scad>;
+include <hollow_out.scad>;
+include <hull_polyline2d.scad>;
+include <hull_polyline3d.scad>;
+include <in_polyline.scad>;
+include <in_shape.scad>;
+include <line2d.scad>;
+include <line3d.scad>;
+include <log.scad>;
+include <m_cumulate.scad>;
+include <m_mirror.scad>;
+include <m_rotation.scad>;
+include <m_scaling.scad>;
+include <m_shearing.scad>;
+include <m_translation.scad>;
+include <midpt_smooth.scad>;
+include <multi_line_text.scad>;
+include <parse_number.scad>;
+include <path_extrude.scad>;
+include <path_scaling_sections.scad>;
+include <paths2sections.scad>;
+include <pie.scad>;
+include <polyline2d.scad>;
+include <polyline3d.scad>;
+include <polysections.scad>;
+include <polytransversals.scad>;
+include <ring_extrude.scad>;
+include <rotate_p.scad>;
+include <rounded_cube.scad>;
+include <rounded_cylinder.scad>;
+include <rounded_extrude.scad>;
+include <rounded_square.scad>;
+include <shape_arc.scad>;
+include <shape_cyclicpolygon.scad>;
+include <shape_ellipse.scad>;
+include <shape_glued2circles.scad>;
+include <shape_path_extend.scad>;
+include <shape_pentagram.scad>;
+include <shape_pie.scad>;
+include <shape_square.scad>;
+include <shape_starburst.scad>;
+include <shape_superformula.scad>;
+include <shape_taiwan.scad>;
+include <shape_trapezium.scad>;
+include <shear.scad>;
+include <sphere_spiral.scad>;
+include <sphere_spiral_extrude.scad>;
+include <split_str.scad>;
+include <starburst.scad>;
+include <stereographic_extrude.scad>;
+include <sub_str.scad>;
+include <torus_knot.scad>;
+include <triangulate.scad>;
+include <trim_shape.scad>;
+include <turtle2d.scad>;
+include <turtle3d.scad>;
+include <voronoi2d.scad>;
+include <voronoi3d.scad>;

docs/lib-along_with.md

@@ -5,9 +5,10 @@ Puts children along the given path. If there's only one child, it will put the c
 ## Parameters
 
 - `points` : The points along the path. 
-- `angles` : Rotate before translate each child. If not given, rotate children automatically according to `points`.
+- `angles` : Rotate before translate each child. If not given, rotate children automatically according to `points` and `method`.
 - `twist` : If given, each child will be twisted before applying each element of `points` and `angles`.
 - `scale` : If given, each child will be scaled before applying each element of `points` and `angles`. It accepts a single value, `[sx, sy]` or `[sx, sy, sz]`.
+- `method` : Which method does `along_with` take to **guess** how to rotate children if `angles` is not specified? It accepts two value, `"AXIS_ANGLE"` (default) and `"EULER_ANGLE"`. See `path_extrude` for more information. **Since:** 1.3.
 
 ## Examples
 
@@ -19,7 +20,7 @@ Puts children along the given path. If there's only one child, it will put the c
 	points = circle_path(radius = 50);
 	
 	along_with(points) 
-	    sphere(5, center = true);
+	    sphere(5);
 
 ![along_with](images/lib-along_with-1.JPG)
 

docs/lib-bend.md

@@ -52,6 +52,4 @@ The arc shape is smoother if the `frags` value is larger.
 
 ![bend](images/lib-bend-3.JPG)
 
-This module is especially useful when you want to create things such as [zentangle bracelet](https://www.thingiverse.com/thing:1569263).
-
-[![zentangle bracelet](http://thingiverse-production-new.s3.amazonaws.com/renders/eb/93/4f/62/1f/3bd1f628e1e566dcb5313035e4f3345b_preview_featured.JPG)](https://www.thingiverse.com/thing:1569263)
+This module is especially useful when you want to create things such as [PNG to pen holder](https://www.thingiverse.com/thing:1589493).

docs/lib-bend_extrude.md

@@ -0,0 +1,36 @@
+# bend_extrude
+
+The purpose of `bend_extrude` is to replace `bend` when you have a 2D shape. `bend_extrude` is faster and doesn't produce jagged edges.
+
+**Since:** 1.3.
+
+## Parameters
+
+- `size` : The size of a square which can contain the target shape.
+- `thickness` : The thinkness used to extrude the shape.
+- `angle` : The central angle of the arc shape. The radius of the arc is calculated automatically.
+- `frags` : Number of fragments. The target shape will be cut into `frags` fragments and recombined into an arc object. The default value is 24.
+
+## Examples
+
+The containing square of the target shape should be laid down on the x-y plane. For example.
+
+	x = 9.25;
+	y = 9.55;
+
+	%square(size = [x, y]);
+	text("A");
+
+![bend_extrude](images/lib-bend_extrude-1.JPG)
+
+Once you have the size of the containing square, you can use it as the `size` argument of the `bend_extrude` module.
+
+	include <bend_extrude.scad>;
+
+	x = 9.25;
+	y = 9.55;
+
+	bend_extrude(size = [x, y], thickness = 1, angle = 270) 
+		text("A");
+
+![bend_extrude](images/lib-bend_extrude-2.JPG)

docs/lib-bijection_offset.md

@@ -0,0 +1,68 @@
+# bijection_offset
+
+Move 2D outlines outward or inward by a given amount. Each point of the offsetted shape is paired with exactly one point of the original shape.
+
+**Since:** 1.2.
+
+## Parameters
+
+- `pts` : Points of a shape.
+- `d` : Amount to offset the shape. When negative, the shape is offset inwards. 
+- `epsilon` : An upper bound on the relative error due to rounding in floating point arithmetic. Default to 0.0001. **Since:** 1.3.
+
+## Examples
+
+	include <bijection_offset.scad>;
+
+	shape = [
+		[15, 0],
+		[15, 30],
+		[0, 20],
+		[-15, 40],
+		[-15, 0]
+	];
+
+	color("red") polygon(bijection_offset(shape, 3));
+	color("orange") polygon(bijection_offset(shape, 2));
+	color("yellow") polygon(bijection_offset(shape, 1));
+	color("green") polygon(shape);
+	color("blue") polygon(bijection_offset(shape, -1));
+	color("indigo") polygon(bijection_offset(shape, -2));
+	color("purple") polygon(bijection_offset(shape, -3));
+
+![bijection_offset](images/lib-bijection_offset-1.JPG)
+
+	include <bijection_offset.scad>;
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <path_extrude.scad>;
+	include <bezier_curve.scad>;
+
+	shape = [
+		[5, 0],
+		[3, 9],
+		[0, 10],    
+		[-5, 0]
+	];
+	offsetted = bijection_offset(shape, 1);
+
+	offsetted2 = bijection_offset(shape, 2);
+	offsetted3 = bijection_offset(shape, 3);
+
+	t_step = 0.05;
+
+	p0 = [0, 0, 0];
+	p1 = [40, 60, 35];
+	p2 = [-50, 70, 0];
+	p3 = [20, 150, -35];
+	p4 = [30, 50, -3];
+
+	path_pts = bezier_curve(t_step, 
+		[p0, p1, p2, p3, p4]
+	);
+
+	path_extrude(concat(offsetted, shape), path_pts, "HOLLOW");
+	path_extrude(concat(offsetted3, offsetted2), path_pts, "HOLLOW");
+
+![bijection_offset](images/lib-bijection_offset-2.JPG)
+

docs/lib-in_polyline.md

@@ -0,0 +1,41 @@
+# in_polyline
+
+Checks wether a point is on a line.
+
+**Since:** 1.3
+
+## Parameters
+
+- `line_pts` : The line points.
+- `pt` : The point to be checked.
+- `epsilon` : An upper bound on the relative error due to rounding in floating point arithmetic. Default to 0.0001.
+
+## Examples
+
+    include <in_polyline.scad>;
+
+    pts = [
+        [0, 0],
+        [10, 0],
+        [10, 10]
+    ];
+
+    echo(in_polyline(pts, [-2, -3]));  // false
+    echo(in_polyline(pts, [5, 0]));    // true
+    echo(in_polyline(pts, [10, 5]));   // true
+    echo(in_polyline(pts, [10, 15]));  // false
+
+----
+
+    include <in_polyline.scad>;
+
+    pts = [
+        [10, 0, 10],
+        [20, 0, 10],
+        [20, 10, 10]
+    ]; 
+
+    echo(in_polyline(pts, [10, 0, 10]));  // true
+    echo(in_polyline(pts, [15, 0, 10]));  // true
+    echo(in_polyline(pts, [15, 1, 10]));  // false
+    echo(in_polyline(pts, [20, 11, 10])); // false    

docs/lib-in_shape.md

@@ -0,0 +1,39 @@
+# in_shape
+
+Checks wether a point is inside a shape.
+
+**Since:** 1.3
+
+## Parameters
+
+- `shapt_pts` : The shape points.
+- `pt` : The point to be checked.
+- `include_edge` : If a point is on the edge of the shape, the function is default to return `false`. If `include_edge` is `true`, the function returns `true`.
+- `epsilon` : An upper bound on the relative error due to rounding in floating point arithmetic. Default to 0.0001.
+
+## Examples
+
+    include <shape_taiwan.scad>;
+    include <in_shape.scad>;
+
+    points = shape_taiwan(30);
+
+    %polygon(points);
+
+    n = 200;
+    xs = rands(-9, 9, n);
+    ys = rands(-16, 16, n);
+
+    pts = [
+        for(i = [0:n - 1]) 
+            let(p = [xs[i], ys[i]]) 
+            if(in_shape(points, p, true))
+            p
+        ]; 
+
+    for(p = pts) {
+        translate(p) 
+            circle(.2); 
+    }
+
+![in_shape](images/lib-in_shape-1.JPG)

docs/lib-m_mirror.md

@@ -0,0 +1,23 @@
+# m_mirror
+
+Generate a 4x4 transformation matrix which can pass into `multmatrix` to mirror the child element on a plane through the origin. 
+
+**Since:** 1.1
+
+## Parameters
+
+- `v` : The normal vector of a plane intersecting the origin through which to mirror the object.
+
+## Examples
+
+	include <m_mirror.scad>;
+
+	rotate([0, 0, 10]) 
+		cube([3, 2, 1]);
+		
+	multmatrix(m_mirror([1, 1, 0]))
+		rotate([0, 0, 10]) 
+			cube([3, 2, 1]);
+
+![m_mirror](images/lib-m_mirror-1.JPG)
+

docs/lib-midpt_smooth.md

@@ -0,0 +1,26 @@
+# midpt_smooth
+
+Given a 2D path, this function constructs a mid-point smoothed version by joining the mid-points of the lines of the path. 
+
+**Since:** 1.3
+
+## Parameters
+
+- `points` : The path points.
+- `n` : Perform mid-point smoothing n times.
+- `closed` : Is the points a 2D shape? If it's `true`, the function takes the last point and the first one to calculate a middle point. Default to `false`.
+
+## Examples
+
+    include <hull_polyline2d.scad>;
+    include <shape_taiwan.scad>;
+    include <bijection_offset.scad>;
+    include <midpt_smooth.scad>;
+
+    taiwan = shape_taiwan(50);  
+    smoothed = midpt_smooth(taiwan, 20, true);
+
+    translate([0, 0, 0]) hull_polyline2d(taiwan, .25); 
+    #translate([10, 0, 0]) hull_polyline2d(smoothed, .25);
+
+![midpt_smooth](images/lib-midpt_smooth-1.JPG)

docs/lib-path_extrude.md

@@ -13,7 +13,8 @@ When using this module, you should use points to represent the 2D shape. If your
 - `triangles` : `"SOLID"` (default), `"HOLLOW"` or user-defined indexes. See example below.
 - `twist` : The number of degrees of through which the shape is extruded.
 - `scale` : Scales the 2D shape by this value over the length of the extrusion. Scale can be a scalar or a vector.
-- `closed` : If the first point and the last point of `path_pts` has the same coordinate, setting `closed` to `true` will connect them automatically.
+- `closed` : If the first point and the last point of `path_pts` has the same coordinate, setting `closed` to `true` will connect them automatically.  
+- `method` : Which method  does `path_extrude` take to **guess** how to generate sections? It accepts two value, `"AXIS_ANGLE"` (default) and `"EULER_ANGLE"`. **Since:** 1.3.
 
 ## Examples
 
@@ -123,6 +124,186 @@ When using this module, you should use points to represent the 2D shape. If your
 
 ![path_extrude](images/lib-path_extrude-3.JPG)
 
+## About `path_extrude` (Important!!)
 
+**`path_extrude` is actually a workaround when you have/provide only path points.**
 
+If you want to extrude a shape along a path precisely, providing enough information about how to rotate sections is necessary. If you want to extrude a shape along a helix, `helix_extrude` is more suitable because it knows how to dig out necessary data for rotating sections precisely.
 
+	include <helix.scad>;
+	include <rotate_p.scad>;
+	include <cross_sections.scad>;
+	include <polysections.scad>;
+	include <helix_extrude.scad>;
+
+	shape_pts = [
+		[0,0],
+		[3, 1],
+		[0, 2]
+	];
+
+	helix_extrude(shape_pts, 
+		radius = 5, 
+		levels = 5, 
+		level_dist = 3,
+		vt_dir = "SPI_UP"
+	);
+
+![path_extrude](images/lib-path_extrude-4.JPG)
+
+If you have only points, what `path_extrude` can do is to **guess** data about rotations. The different algorithm will dig out different data. For example:
+
+	include <helix.scad>;
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <helix.scad>;
+	include <path_extrude.scad>;
+
+	shape_pts = [
+		[0,0],
+		[3, 1],
+		[0, 2]
+	];
+
+	points = helix( 
+		radius = 5, 
+		levels = 5, 
+		level_dist = 3,
+		vt_dir = "SPI_UP"
+	);
+
+	path_extrude(shape_pts, points);
+
+![path_extrude](images/lib-path_extrude-5.JPG)
+
+You might think this is wrong. Actually, it's not. It's the correct/default behavior of `path_extrude`. Because **you don't provide other information**, what `path_extrude` can do is to **guess** how to generate sections from points. You think it's a bug in `path_extrude` because your brain has information that path points do not provide.
+
+ The `method` parameter is default to `"AXIS_ANGLE"`, a way to guess information from points. It accepts `"EULER_ANGLE"`, too. 
+
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <helix.scad>;
+	include <path_extrude.scad>;
+
+	shape_pts = [
+		[0,0],
+		[3, 1],
+		[0, 2]
+	];
+
+	points = helix( 
+		radius = 5, 
+		levels = 5, 
+		level_dist = 3,
+		vt_dir = "SPI_UP"
+	);
+
+	path_extrude(shape_pts, points, method = "EULER_ANGLE");
+
+![path_extrude](images/lib-path_extrude-6.JPG)
+
+`"EULER_ANGLE"` generates the same section at the same point. You might think the model is correct. But, that's because what it guesses from points just match your expectation. 
+
+`"EULER_ANGLE"` will generate an abrupt when the path is exactly vertical. [The problem happened in (older) Blender, too.](https://download.blender.org/documentation/htmlI/ch09s04.html) 
+
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <path_extrude.scad>;
+
+	shape_pts = [[5, -5], [5, 5], [-5, 5], [-5, -5]];
+
+	path_pts = [
+		[20, 20, 0], 
+		[18.2, 18.2, 2], 
+		[16.8, 16.8, 4], 
+		[15.8, 15.8, 6], 
+		[15.2, 15.2, 8], 
+		[15, 15, 10], 
+		[15.2, 15.2, 12], 
+		[15.8, 15.8, 14], 
+		[16.8, 16.8, 16], 
+		[18.2, 18.2, 18], 
+		[20, 20, 20]
+	];
+
+	path_extrude(shape_pts, path_pts, method = "EULER_ANGLE");
+
+![path_extrude](images/lib-path_extrude-7.JPG)
+
+The problem doesn't happen when `method` is `"AXIS_ANGLE"`.
+
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <path_extrude.scad>;
+
+	shape_pts = [[5, -5], [5, 5], [-5, 5], [-5, -5]];
+
+	path_pts = [
+		[20, 20, 0], 
+		[18.2, 18.2, 2], 
+		[16.8, 16.8, 4], 
+		[15.8, 15.8, 6], 
+		[15.2, 15.2, 8], 
+		[15, 15, 10], 
+		[15.2, 15.2, 12], 
+		[15.8, 15.8, 14], 
+		[16.8, 16.8, 16], 
+		[18.2, 18.2, 18], 
+		[20, 20, 20]
+	];
+
+	path_extrude(shape_pts, path_pts, method = "AXIS_ANGLE");
+
+![path_extrude](images/lib-path_extrude-8.JPG)
+
+So, which is the correct method? Both methods are correct when you provide only points. `method` is just a way you tell `path_extrude` how to guess more information when extruding. 
+
+`"EULER_ANGLE"` will generate an abrupt when the path is exactly vertical. Some users might think it's a bug so `"AXIS_ANGLE"` is the default value. 
+
+`"EULER_ANGLE"`, however, generates the same section at the same point. This means that you don't have to adjust sections if you want to extrude along a closed path. It's an advantage when extruding. For example:
+
+	include <shape_pentagram.scad>;
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <path_extrude.scad>;
+	include <torus_knot.scad>;
+
+	p = 2;
+	q = 3;
+	phi_step = 0.05;
+	star_radius = 0.5;
+
+	pts = torus_knot(p, q, phi_step);
+
+	shape_pentagram_pts = shape_pentagram(star_radius);
+
+	// not closed perfectly
+	translate([-8, 0, 0]) path_extrude(
+		shape_pentagram_pts, 
+		concat(pts, [pts[0]]), 
+		closed = true,
+		method = "AXIS_ANGLE"
+	);
+
+	// adjust it 
+	path_extrude(
+		shape_pentagram_pts, 
+		concat(pts, [pts[0]]), 
+		closed = true,
+		twist = 188,
+		method = "AXIS_ANGLE"
+	);
+
+	// "EULER_ANGLE" is easy in this situation
+	translate([0, 8, 0]) path_extrude(
+		shape_pentagram_pts, 
+		concat(pts, [pts[0]]), 
+		closed = true,
+		method = "EULER_ANGLE"
+	);
+
+![path_extrude](images/lib-path_extrude-9.JPG)	
+
+Both methods are useful. If `"AXIS_ANGLE"` doesn't guess out what you want, choose `"EULER_ANGLE"`, and vice versa.
+
+For more information, see [#issue 3](https://github.com/JustinSDK/dotSCAD/issues/3) and [#issue 5](https://github.com/JustinSDK/dotSCAD/issues/5).

docs/lib-path_scaling_sections.md

@@ -0,0 +1,133 @@
+# path_scaling_sections
+
+Given an edge path with the first point at the outline of a shape. This function uses the path to calculate scaling factors and returns all scaled sections in the reversed order of the edge path. Combined with the `polysections` module, you can create an extrusion with the path as an edge. 
+
+In order to control scaling factors easily, I suggest using `[x, 0, 0]` as the first point and keeping y = 0 while building the edge path.
+
+You can use any point as the first point of the edge path. Just remember that your edge path radiates from the origin.
+
+**Since:** 1.2.
+
+## Parameters
+
+- `shape_pts` : A list of points represent a shape.
+- `edge_path` : A list of points represent the edge path.
+
+## Examples
+
+	include <hull_polyline3d.scad>;
+	include <shape_taiwan.scad>;
+	include <path_scaling_sections.scad>;
+	include <m_scaling.scad>;
+	include <polysections.scad>;
+
+	taiwan = shape_taiwan(100);
+	fst_pt = [13, 0, 0];
+
+	edge_path = [
+		fst_pt,
+		fst_pt + [0, 0, 10],
+		fst_pt + [10, 0, 20],
+		fst_pt + [8, 0, 30],
+		fst_pt + [12, 0, 40],
+		fst_pt + [0, 0, 50],
+		fst_pt + [0, 0, 60]
+	];
+
+	#hull_polyline3d(edge_path);
+	polysections(path_scaling_sections(taiwan, edge_path));
+
+![path_scaling_sections](images/lib-path_scaling_sections-1.JPG)
+
+	include <hull_polyline3d.scad>;
+	include <shape_taiwan.scad>;
+	include <path_scaling_sections.scad>;
+	include <m_scaling.scad>;
+	include <polysections.scad>;
+	include <bezier_curve.scad>;
+
+
+	taiwan = shape_taiwan(100);
+	fst_pt = [13, 0, 0];
+
+	edge_path = bezier_curve(0.05, [
+		fst_pt,
+		fst_pt + [0, 0, 10],
+		fst_pt + [10, 0, 20],
+		fst_pt + [8, 0, 30],
+		fst_pt + [12, 0, 40],
+		fst_pt + [0, 0, 50],
+		fst_pt + [0, 0, 60]
+	]);
+
+	#hull_polyline3d(edge_path);
+	polysections(path_scaling_sections(taiwan, edge_path));
+
+![path_scaling_sections](images/lib-path_scaling_sections-2.JPG)
+
+	include <shape_taiwan.scad>;
+	include <path_scaling_sections.scad>;
+	include <m_scaling.scad>;
+	include <polysections.scad>;
+	include <bezier_curve.scad>;
+	include <rotate_p.scad>;
+
+	taiwan = shape_taiwan(100);
+	fst_pt = [13, 0, 0];
+
+	edge_path = bezier_curve(0.05, [
+		fst_pt,
+		fst_pt + [0, 0, 10],
+		fst_pt + [10, 0, 20],
+		fst_pt + [8, 0, 30],
+		fst_pt + [12, 0, 40],
+		fst_pt + [0, 0, 50],
+		fst_pt + [0, 0, 60]
+	]);
+
+	leng = len(edge_path);
+	twist = -90;
+	twist_step = twist / leng;
+	sections = path_scaling_sections(taiwan, edge_path);
+
+	rotated_sections = [
+		for(i = [0:leng - 1]) 
+		[
+			for(p = sections[i]) 
+				rotate_p(p, twist_step * i)        
+		]
+	];
+
+	polysections(rotated_sections);
+
+![path_scaling_sections](images/lib-path_scaling_sections-3.JPG)	
+
+	include <hull_polyline3d.scad>;
+	include <shape_taiwan.scad>;
+	include <path_scaling_sections.scad>;
+	include <m_scaling.scad>;
+	include <polysections.scad>;
+    include <rotate_p.scad>;
+
+	taiwan = shape_taiwan(100);
+
+    /* 
+	    You can use any point as the first point of the edge path.
+		Just remember that your edge path radiates from the origin.
+	*/
+	fst_pt = [taiwan[0][0], taiwan[0][1], 0];//[13, 0, 0];
+    a = atan2(fst_pt[1], fst_pt[0]);
+	edge_path = [
+		fst_pt,
+		fst_pt + rotate_p([0, 0, 10], a),
+		fst_pt + rotate_p([10, 0, 20], a),
+		fst_pt + rotate_p([8, 0, 30], a),
+		fst_pt + rotate_p([10, 0, 40], a),
+		fst_pt + rotate_p([0, 0, 50], a),
+		fst_pt + rotate_p([0, 0, 60], a)
+	];
+
+	#hull_polyline3d(edge_path);
+	polysections(path_scaling_sections(taiwan, edge_path));
+
+![path_scaling_sections](images/lib-path_scaling_sections-4.JPG)

docs/lib-shape_starburst.md

@@ -6,7 +6,7 @@ Returns shape points of a star. They can be used with xxx_extrude modules of dot
 
 - `r1` : The outer radius of the starburst. 
 - `r2` : The inner radius of the starburst.
-- `n`  : The number of vertices. 
+- `n`  : The burst number. 
 
 
 ## Examples

docs/lib-shape_taiwan.md

@@ -5,6 +5,7 @@ Returns shape points of [Taiwan](https://www.google.com.tw/maps?q=taiwan&um=1&ie
 ## Parameters
 
 - `h` : The height of Taiwan.
+- `distance` : Used for simplifying the shape. If the distance between a point and its previous points is not greater than `distance`, the point will be kept. Default to 0. **Since:** 1.3.
 
 ## Examples
 

docs/lib-starburst.md

@@ -0,0 +1,24 @@
+# starburst 
+
+A 3D version of `shape_starburst`. 
+
+**Since:** 1.2.
+
+## Parameters
+
+- `r1` : The outer radius of the starburst. 
+- `r2` : The inner radius of the starburst.
+- `n`  : The number of vertices. 
+- `height` : The height of the starburst.
+
+## Examples
+
+    include <starburst.scad>;
+
+	starburst(10, 5, 5, 5);
+	translate([20, 0, 0]) starburst(10, 5, 6, 5);
+	translate([40, 0, 0]) starburst(10, 5, 12, 10);
+	translate([60, 0, 0]) starburst(10, 5, 4, 3);
+
+![starburst](images/lib-starburst-1.JPG)
+

docs/lib-torus_knot.md

@@ -0,0 +1,40 @@
+# torus_knot 
+
+Generate a path of [The (p,q)-torus knot](https://en.wikipedia.org/wiki/Torus_knot).
+
+**Since:** 1.2.
+
+![torus_knot](images/lib-torus_knot-1.JPG)
+
+## Parameters
+
+- `p` : The p parameter of The (p,q)-torus knot. 
+- `q` : The q parameter of The (p,q)-torus knot. 
+- `phi_step` : The amount when increasing phi. 
+
+## Examples
+
+	include <shape_pentagram.scad>;
+	include <rotate_p.scad>;
+	include <polysections.scad>;
+	include <path_extrude.scad>;
+	include <torus_knot.scad>;
+
+	p = 2;
+	q = 3;
+	phi_step = 0.05;
+	star_radius = 0.5;
+
+	pts = torus_knot(p, q, phi_step);
+
+	shape_pentagram_pts = shape_pentagram(star_radius);
+
+	path_extrude(
+		shape_pentagram_pts, 
+		concat(pts, [pts[0]]), 
+		closed = true,
+		method = "EULER_ANGLE"
+	);
+
+![torus_knot](images/lib-torus_knot-2.JPG)
+

docs/lib-triangulate.md

@@ -0,0 +1,42 @@
+# triangulate
+
+Given a 2D shape. This function performs a simple polygon triangulation algorithm and returns the indices of each triangle. 
+
+**Since:** 1.3.
+
+## Parameters
+
+- `shape_pts` : The shape points.
+- `epsilon` : An upper bound on the relative error due to rounding in floating point arithmetic. Default to 0.0001.
+
+## Examples
+
+    include <triangulate.scad>; 
+
+    shape = [
+        [0, 0],
+        [10, 0],
+        [12, 5],
+        [5, 10],
+        [10, 15],
+        [0, 20],
+        [-5, 18],
+        [-18, 3],
+        [-4, 10]
+    ];
+
+    tris = triangulate(shape);
+
+    difference() {
+        polygon(shape);
+
+        for(tri = tris) {
+            offset(-.2) 
+                polygon([for(idx = tri) shape[idx]]);
+            
+        }
+    }
+
+
+![triangulate](images/lib-triangulate-1.JPG)
+

docs/lib-trim_shape.md

@@ -0,0 +1,33 @@
+# trim_shape
+
+Given a tangled-edge shape. This function trims the shape to a non-tangled shape. It's intended to be a helper function after using `bijection_offset`. 
+
+**Since:** 1.3.
+
+## Parameters
+
+- `shape_pts` : The shape points.
+- `from` : The index of the start point you want to trim.
+- `to` : The index of the last point you want to trim.
+- `epsilon` : An upper bound on the relative error due to rounding in floating point arithmetic. Default to 0.0001.
+
+## Examples
+
+    include <hull_polyline2d.scad>;
+    include <trim_shape.scad>;
+    include <shape_taiwan.scad>;
+    include <bijection_offset.scad>;
+    include <midpt_smooth.scad>;
+
+    taiwan = shape_taiwan(50);
+    offseted = bijection_offset(taiwan, -2);
+    trimmed = trim_shape(offseted, 3, len(offseted) - 6);
+    smoothed = midpt_smooth(trimmed, 3);
+
+    #hull_polyline2d(taiwan, .1); 
+    %translate([25, 0, 0]) 
+        hull_polyline2d(offseted, .2);
+    hull_polyline2d(smoothed, .1); 
+
+![trim_shape](images/lib-trim_shape-1.JPG)
+

docs/lib-voronoi2d.md

@@ -0,0 +1,47 @@
+# voronoi2d
+
+Creats a [Voronoi diagram](https://en.wikipedia.org/wiki/Voronoi_diagram). The initial region for each cell is calculated automatically from the given points by the following code: 
+
+    xs = [for(p = points) p[0]];
+    ys = [for(p = points) abs(p[1])];
+    region_size = max([(max(xs) -  min(xs) / 2), (max(ys) -  min(ys)) / 2]);    
+
+**Since:** 1.3.
+
+## Parameters
+
+- `points` : Points for each cell. 
+- `spacing` : Distance between cells. Default to 1.
+- `r`, `delta`, `chamfer` : The outlines of each cell can be moved outward or inward. These parameters have the same effect as [`offset`](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Transformations#offset). 
+- `region_type` : The initial shape for each cell can be `"square"` or `"circle"`. Default to `"square"`.
+
+## Examples
+
+    include <voronoi2d.scad>;
+
+    xs = rands(-20, 20, 50);
+    ys = rands(-20, 20, 50);
+
+    points = [for(i = [0:len(xs) - 1]) [xs[i], ys[i]]];
+
+    voronoi2d(points);
+    translate([60, 0, 0]) 
+        voronoi(points, region_type = "circle");
+
+![voronoi2d](images/lib-voronoi2d-1.JPG)
+
+    include <voronoi2d.scad>;
+    include <hollow_out.scad>;
+
+    xs = rands(0, 40, 50);
+    ys = rands(0, 20, 50);
+
+    points = [for(i = [0:len(xs) - 1]) [xs[i], ys[i]]];
+
+    difference() {
+        square([40, 20]);
+        voronoi2d(points);
+    }
+    hollow_out(shell_thickness = 1) square([40, 20]);
+    
+![voronoi2d](images/lib-voronoi2d-2.JPG)

docs/lib-voronoi3d.md

@@ -0,0 +1,60 @@
+# voronoi3d
+
+Creats a 3D version of [Voronoi diagram](https://en.wikipedia.org/wiki/Voronoi_diagram). The initial space for each cell is calculated automatically from the given points by the following code: 
+
+    xs = [for(p = points) p[0]];
+    ys = [for(p = points) abs(p[1])];
+    zs = [for(p = points) abs(p[2])];
+    space_size = max([(max(xs) -  min(xs) / 2), (max(ys) -  min(ys)) / 2, (max(zs) -  min(zs)) / 2]);
+    // cube([space_size, space_size * 2, space_size * 2]);
+
+The preview or rendering of 3D Voronoi is slow. If you want to use this module, render and export the 3D Voronoi model first. Then, `import` the model to do what you want.
+
+**Since:** 1.3.
+
+## Parameters
+
+- `points` : Points for each cell. 
+- `spacing` : Distance between cells. Default to 1.
+
+## Examples
+
+    include <voronoi3d.scad>;
+
+    r = 30;
+
+    zas = rands(0, 359, 12);
+    yas = rands(0, 179, 12);
+
+    points = [
+        for(i = [0:len(zas) - 1])
+        [
+            r * cos(yas[i]) * cos(zas[i]), 
+            r * cos(yas[i]) * sin(zas[i]), 
+            r * sin(yas[i])
+        ]
+    ];
+
+    #for(pt = points) {
+        translate(pt) cube(1);
+    }
+
+    intersection() {
+        sphere(r);
+        voronoi3d(points);
+    }
+
+![voronoi3d](images/lib-voronoi3d-1.JPG)
+
+If you render, export and save the previous model as `voronoi3d.stl`, the following code will generate a Voronoi sphere.
+
+    r = 30;
+    thickness = 2;
+
+    difference() {
+        sphere(r);
+        scale(1.01) import("voronoi3d.stl");
+        sphere(r - thickness);
+    }
+    
+![voronoi3d](images/lib-voronoi3d-2.JPG)

src/__private__/__half_trapezium.scad

@@ -44,9 +44,9 @@ function __br_corner(frags, b_ang, l1, l2, h, round_r) =
 
 function __half_trapezium(length, h, round_r) =
     let(
-        is_vt = __is_vector(length),
-        l1 = is_vt ? length[0] : length,
-        l2 = is_vt ? length[1] : length,
+        is_flt = __is_float(length),
+        l1 = is_flt ? length : length[0],
+        l2 = is_flt ? length : length[1],
         frags = __frags(round_r),
         b_ang = atan2(h, l1 - l2),
         br_corner = __br_corner(frags, b_ang, l1, l2, h, round_r),

src/__private__/__in_line.scad

@@ -0,0 +1,8 @@
+function __in_line(line_pts, pt, epsilon = 0.0001) =
+    let(
+        pts = len(line_pts[0]) == 2 ? [for(p = line_pts) __to3d(p)] : line_pts,
+        pt3d = len(pt) == 2 ? __to3d(pt) : pt,
+        v1 = pts[0] - pt3d, 
+        v2 = pts[1] - pt3d
+    )
+    (norm(cross(v1, v2)) < epsilon) && ((v1 * v2) <= epsilon);

src/__private__/__line_intersection.scad

@@ -0,0 +1,12 @@
+function __line_intersection(line_pts1, line_pts2, epsilon = 0.0001) = 
+    let(
+        a1 = line_pts1[0],
+        a2 = line_pts1[1],
+        b1 = line_pts2[0],
+        b2 = line_pts2[1],
+        a = a2 - a1, 
+        b = b2 - b1, 
+        s = b1 - a1
+    )
+    abs(cross(a, b)) < epsilon ? [] :  // they are parallel or conincident edges
+        a1 + a * cross(s, b) / cross(a, b);

src/__private__/__lines_from.scad

@@ -0,0 +1,7 @@
+function __lines_from(pts, closed = false) = 
+    let(leng = len(pts))
+    concat(
+        [for(i = [0:leng - 2]) [pts[i], pts[i + 1]]], 
+        closed ? [[pts[len(pts) - 1], pts[0]]] : []
+    );
+    

src/__private__/__shape_arc.scad

@@ -6,7 +6,7 @@ function __shape_arc(radius, angle, width, width_mode = "LINE_CROSS") =
         frags = __frags(radius),
         a_step = 360 / frags,
         half_a_step = a_step / 2,
-        angles = __is_vector(angle) ? angle : [0, angle],
+        angles = __is_float(angle) ? [0, angle] : angle,
         m = floor(angles[0] / a_step) + 1,
         n = floor(angles[1] / a_step),
         r_outer = radius + w_offset[0],

src/__private__/__shape_pie.scad

@@ -3,7 +3,7 @@ function __shape_pie(radius, angle) =
         frags = __frags(radius),
         a_step = 360 / frags,
         leng = radius * cos(a_step / 2),
-        angles = __is_vector(angle) ? angle : [0:angle],
+        angles = __is_float(angle) ? [0:angle] : angle,
         m = floor(angles[0] / a_step) + 1,
         n = floor(angles[1] / a_step),
         edge_r_begin = leng / cos((m - 0.5) * a_step - angles[0]),

src/__private__/__to_ang_vect.scad

@@ -0,0 +1,7 @@
+function __to_3_elems_ang_vect(a) =
+     let(leng = len(a))
+     leng == 3 ? a : (
+         leng == 2 ? [a[0], a[1], 0] :  [a[0], 0, 0]
+     );
+
+function __to_ang_vect(a) = __is_float(a) ? [0, 0, a] : __to_3_elems_ang_vect(a);

src/along_with.scad

@@ -9,31 +9,42 @@
 **/ 
  
 include <__private__/__angy_angz.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__to3d.scad>;
 
 // Becuase of improving the performance, this module requires m_rotation.scad which doesn't require in dotSCAD 1.0. 
 // For backward compatibility, I directly include m_rotation here.
 include <m_rotation.scad>;
 
-module along_with(points, angles, twist = 0, scale = 1.0) {
+module along_with(points, angles, twist = 0, scale = 1.0, method = "AXIS_ANGLE") {
     leng_points = len(points);
     leng_points_minus_one = leng_points - 1;
     twist_step_a = twist / leng_points;
-        
-    function scale_step() =
-        let(s =  (scale - 1) / leng_points_minus_one)
-        [s, s, s];
 
-    scale_step_vt = __is_vector(scale) ? 
+    angles_defined = angles != undef;
+
+    scale_step_vt = __is_float(scale) ? 
+        scale_step() :
         [
             (scale[0] - 1) / leng_points_minus_one, 
             (scale[1] - 1) / leng_points_minus_one,
             scale[2] == undef ? 0 : (scale[2] - 1) / leng_points_minus_one
-        ] : scale_step(); 
+        ]; 
+
+
+    function scale_step() =
+        let(s =  (scale - 1) / leng_points_minus_one)
+        [s, s, s];
+
+
+    /* 
+         Sadly, children(n) cannot be used with inner modules 
+         so I have to do things in the first level. Ugly!!
+    */
+
+    // >>> begin: modules and functions for "AXIS-ANGLE"
 
     // get rotation matrice for sections
-
     identity_matrix = [
         [1, 0, 0, 0],
         [0, 1, 0, 0],
@@ -41,19 +52,19 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
         [0, 0, 0, 1]
     ];    
 
-    function local_ang_vects(j) = 
-        j == 0 ? [] : local_ang_vects_sub(j);
+    function axis_angle_local_ang_vects(j) = 
+        j == 0 ? [] : axis_angle_local_ang_vects_sub(j);
     
-    function local_ang_vects_sub(j) =
+    function axis_angle_local_ang_vects_sub(j) =
         let(
             vt0 = points[j] - points[j - 1],
             vt1 = points[j + 1] - points[j],
             a = acos((vt0 * vt1) / (norm(vt0) * norm(vt1))),
             v = cross(vt0, vt1)
         )
-        concat([[a, v]], local_ang_vects(j - 1));
+        concat([[a, v]], axis_angle_local_ang_vects(j - 1));
 
-    function cumulated_rot_matrice(i, rot_matrice) = 
+    function axis_angle_cumulated_rot_matrice(i, rot_matrice) = 
         let(
             leng_rot_matrice = len(rot_matrice),
             leng_rot_matrice_minus_one = leng_rot_matrice - 1,
@@ -66,13 +77,13 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
                    rot_matrice[leng_rot_matrice_minus_one], 
                    rot_matrice[leng_rot_matrice_minus_two] * rot_matrice[leng_rot_matrice_minus_one]
                ] 
-               : cumulated_rot_matrice_sub(i, rot_matrice)
+               : axis_angle_cumulated_rot_matrice_sub(i, rot_matrice)
             )
         );
 
-    function cumulated_rot_matrice_sub(i, rot_matrice) = 
+    function axis_angle_cumulated_rot_matrice_sub(i, rot_matrice) = 
         let(
-            matrice = cumulated_rot_matrice(i + 1, rot_matrice),
+            matrice = axis_angle_cumulated_rot_matrice(i + 1, rot_matrice),
             curr_matrix = rot_matrice[i],
             prev_matrix = matrice[len(matrice) - 1]
         )
@@ -80,7 +91,7 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
 
     // align modules
 
-    module align_with_pts_angles(i) {
+    module axis_angle_align_with_pts_angles(i) {
         translate(points[i]) 
             rotate(angles[i])
                 rotate(twist_step_a * i) 
@@ -88,8 +99,9 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
                             children(0);
     }
 
-    module align_with_pts_init(a, s) {
-        angleyz = __angy_angz(points[0], points[1]);
+    module axis_angle_align_with_pts_init(a, s) {
+        angleyz = __angy_angz(__to3d(points[0]), __to3d(points[1]));
+
         rotate([0, -angleyz[0], angleyz[1]])
             rotate([90, 0, -90])
                 rotate(a)
@@ -97,51 +109,102 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
                         children(0);
     }
     
-    module align_with_pts_local_rotate(j, init_a, init_s, cumu_rot_matrice) {
+    module axis_angle_align_with_pts_local_rotate(j, init_a, init_s, cumu_rot_matrice) {
         if(j == 0) {  // first child
-            align_with_pts_init(init_a, init_s) 
+            axis_angle_align_with_pts_init(init_a, init_s) 
                 children(0);
         }
         else {
             multmatrix(cumu_rot_matrice[j - 1])
-                align_with_pts_init(init_a, init_s) 
+                axis_angle_align_with_pts_init(init_a, init_s) 
                     children(0);
         }
     } 
 
-    if(angles != undef) {
+    // <<< end: modules and functions for "AXIS-ANGLE"
+
+
+    // >>> begin: modules and functions for "EULER-ANGLE"
+
+    function _euler_angle_path_angles(pts, end_i, i = 0) = 
+        i == end_i ?
+                [] : 
+                concat(
+                    [__angy_angz(pts[i], pts[i + 1])], 
+                    _euler_angle_path_angles(pts, end_i, i + 1)
+                );
+            
+    function euler_angle_path_angles(children) = 
+       let(
+           pts = len(points[0]) == 3 ? points : [for(pt = points) __to3d(pt)],
+           end_i = children == 1 ? leng_points_minus_one : children - 1,
+           angs = _euler_angle_path_angles(pts, end_i)
+        )
+       concat(
+           [[0, -angs[0][0], angs[0][1]]], 
+           [for(a = angs) [0, -a[0], a[1]]]
+       );
+
+    module euler_angle_align(i, angs) {
+        translate(points[i]) 
+            rotate(angs[i])
+                rotate(angles_defined ? [0, 0, 0] : [90, 0, -90])
+                    rotate(twist_step_a * i) 
+                         scale([1, 1, 1] + scale_step_vt * i) 
+                             children(0);
+    }
+
+    // <<< end: modules and functions for "EULER-ANGLE"
+
+    if(method == "AXIS_ANGLE") {
+        if(angles_defined) {
+            if($children == 1) { 
+                for(i = [0:leng_points_minus_one]) {
+                    axis_angle_align_with_pts_angles(i) children(0);
+                }
+            } else {
+                for(i = [0:min(leng_points, $children) - 1]) {
+                    axis_angle_align_with_pts_angles(i) children(i);
+                }
+            }
+        }
+        else {
+            cumu_rot_matrice = axis_angle_cumulated_rot_matrice(0, [
+                for(ang_vect = axis_angle_local_ang_vects(leng_points - 2)) 
+                    m_rotation(ang_vect[0], ang_vect[1])
+            ]);
+
+            translate(points[0])
+                axis_angle_align_with_pts_local_rotate(0, 0, [1, 1, 1], cumu_rot_matrice)
+                    children(0); 
+
+            if($children == 1) { 
+                for(i = [0:leng_points - 2]) {
+                    translate(points[i + 1])
+                        axis_angle_align_with_pts_local_rotate(i, i * twist_step_a, [1, 1, 1] + scale_step_vt * i, cumu_rot_matrice)
+                            children(0);          
+                }          
+            } else {
+                for(i = [0:min(leng_points, $children) - 2]) {
+                    translate(points[i + 1])
+                        axis_angle_align_with_pts_local_rotate(i, i * twist_step_a, [1, 1, 1] + scale_step_vt * i, cumu_rot_matrice)
+                            children(i + 1);   
+                }
+            }
+        }
+    }
+    else if(method == "EULER_ANGLE") {
+        angs = angles_defined ? angles : euler_angle_path_angles($children);
+
         if($children == 1) { 
             for(i = [0:leng_points_minus_one]) {
-                align_with_pts_angles(i) children(0);
+                euler_angle_align(i, angs) children(0);
             }
+            
         } else {
             for(i = [0:min(leng_points, $children) - 1]) {
-                align_with_pts_angles(i) children(i);
+                euler_angle_align(i, angs) children(i);
             }
-        }
-    }
-    else {
-        cumu_rot_matrice = cumulated_rot_matrice(0, [
-            for(ang_vect = local_ang_vects(leng_points - 2)) 
-                m_rotation(ang_vect[0], ang_vect[1])
-        ]);
-
-        translate(points[0])
-            align_with_pts_local_rotate(0, 0, [1, 1, 1], cumu_rot_matrice)
-                children(0); 
-
-        if($children == 1) { 
-            for(i = [0:leng_points - 2]) {
-                translate(points[i + 1])
-                    align_with_pts_local_rotate(i, i * twist_step_a, [1, 1, 1] + scale_step_vt * i, cumu_rot_matrice)
-                        children(0);          
-            }          
-        } else {
-            for(i = [0:min(leng_points, $children) - 2]) {
-                translate(points[i + 1])
-                    align_with_pts_local_rotate(i, i * twist_step_a, [1, 1, 1] + scale_step_vt * i, cumu_rot_matrice)
-                        children(i + 1);   
-            }
-        }
+        }         
     }
 }

src/arc.scad

@@ -9,7 +9,7 @@
 **/ 
 
 include <__private__/__frags.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__ra_to_xy.scad>;
 include <__private__/__edge_r.scad>;
 include <__private__/__shape_arc.scad>;

src/arc_path.scad

@@ -9,7 +9,7 @@
 **/ 
 
 include <__private__/__frags.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__ra_to_xy.scad>;
 include <__private__/__edge_r.scad>;
 
@@ -17,7 +17,7 @@ function arc_path(radius, angle) =
     let(
         frags = __frags(radius),
         a_step = 360 / frags,
-        angles = __is_vector(angle) ? angle : [0, angle],
+        angles = __is_float(angle) ? [0, angle] : angle,
         m = floor(angles[0] / a_step) + 1,
         n = floor(angles[1] / a_step),
         points = concat([__ra_to_xy(__edge_r_begin(radius, angles[0], a_step, m), angles[0])],

src/bend.scad

@@ -43,7 +43,7 @@ module bend(size, angle, frags = 24) {
             }
     }
 
-    for(i = [0 : frags - 1]) {
+    rotate(90) for(i = [0 : frags - 1]) {
         rotate(i * frag_angle + half_frag_angle) 
             get_frag(i) 
                 children();  

src/bend_extrude.scad

@@ -0,0 +1,41 @@
+/**
+* bend_extrude.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-bend_extrude.html
+*
+**/
+
+module bend_extrude(size, thickness, angle, frags = 24) {
+    x = size[0];
+    y = size[1];
+    frag_width = x / frags ;
+    frag_angle = angle / frags;
+    half_frag_width = 0.5 * frag_width;
+    half_frag_angle = 0.5 * frag_angle;
+    r = half_frag_width / sin(half_frag_angle);
+    s =  (r - thickness) / r;
+    
+    module get_frag(i) {
+        offsetX = i * frag_width;
+        linear_extrude(thickness, scale = [s, 1]) 
+            translate([-offsetX - half_frag_width, 0, 0]) 
+                intersection() {
+                    translate([x, 0, 0]) mirror([1, 0, 0]) children();
+                    translate([offsetX, 0, 0]) 
+                        square([frag_width, y]);
+                }
+    }
+
+    offsetY = -r * cos(half_frag_angle) ;
+    rotate([180, 0, 180]) for(i = [0 : frags - 1]) {
+       rotate(i * frag_angle + half_frag_angle) 
+            translate([0, offsetY, 0])
+                rotate([-90, 0, 0]) 
+                    get_frag(i) 
+                        children();  
+    }
+
+}

src/bijection_offset.scad

@@ -0,0 +1,65 @@
+/**
+* bijection_offset.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-bijection_offset.html
+*
+**/
+
+include <__private__/__lines_from.scad>;
+include <__private__/__line_intersection.scad>;
+    
+function _bijection_inward_edge_normal(edge) =  
+    let(
+        pt1 = edge[0],
+        pt2 = edge[1],
+        dx = pt2[0] - pt1[0],
+        dy = pt2[1] - pt1[1],
+        edge_leng = norm([dx, dy])
+    )
+    [-dy / edge_leng, dx / edge_leng];
+
+function _bijection_outward_edge_normal(edge) = -1 * _bijection_inward_edge_normal(edge);
+
+function _bijection_offset_edge(edge, dx, dy) = 
+    let(
+        pt1 = edge[0],
+        pt2 = edge[1],
+        dxy = [dx, dy]
+    )
+    [pt1 + dxy, pt2 + dxy];
+    
+function _bijection__bijection_offset_edges(edges, d) = 
+    [ 
+        for(edge = edges)
+        let(
+            ow_normal = _bijection_outward_edge_normal(edge),
+            dx = ow_normal[0] * d,
+            dy = ow_normal[1] * d
+        )
+        _bijection_offset_edge(edge, dx, dy)
+    ];
+
+function bijection_offset(pts, d, epsilon = 0.0001) = 
+    let(
+        es = __lines_from(pts, true), 
+        offset_es = _bijection__bijection_offset_edges(es, d),
+        leng = len(offset_es),
+        last_p = __line_intersection(offset_es[leng - 1], offset_es[0], epsilon)
+    )
+    concat(
+        [
+            for(i = [0:leng - 2]) 
+            let(
+                this_edge = offset_es[i],
+                next_edge = offset_es[i + 1],
+                p = __line_intersection(this_edge, next_edge, epsilon)
+            )
+            // p == p to avoid [nan, nan], because [nan, nan] != [nan, nan]
+            if(p != [] && p == p) p
+        ],
+        last_p != [] && last_p == last_p ? [last_p] : []
+    );
+    

src/cross_sections.scad

@@ -9,16 +9,17 @@
 **/
 
 include <__private__/__to3d.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 
 function cross_sections(shape_pts, path_pts, angles, twist = 0, scale = 1.0) =
     let(
         len_path_pts_minus_one = len(path_pts) - 1,
         sh_pts = len(shape_pts[0]) == 3 ? shape_pts : [for(p = shape_pts) __to3d(p)],
         pth_pts = len(path_pts[0]) == 3 ? path_pts : [for(p = path_pts) __to3d(p)],
-        scale_step_vt = __is_vector(scale) ? 
-            [(scale[0] - 1) / len_path_pts_minus_one, (scale[1] - 1) / len_path_pts_minus_one] :
-            [(scale - 1) / len_path_pts_minus_one, (scale - 1) / len_path_pts_minus_one],
+        scale_step_vt = __is_float(scale) ? 
+            [(scale - 1) / len_path_pts_minus_one, (scale - 1) / len_path_pts_minus_one] :
+            [(scale[0] - 1) / len_path_pts_minus_one, (scale[1] - 1) / len_path_pts_minus_one]
+            ,
         scale_step_x = scale_step_vt[0],
         scale_step_y = scale_step_vt[1],
         twist_step = twist / len_path_pts_minus_one

src/crystal_ball.scad

@@ -9,10 +9,10 @@
 **/ 
 
 include <__private__/__nearest_multiple_of_4.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 
 module crystal_ball(radius, theta = 360, phi = 180) {
-    phis = __is_vector(phi) ? phi : [0, phi];
+    phis = __is_float(phi) ? [0, phi] : phi;
     
     frags = __frags(radius);
 

src/helix.scad

@@ -8,14 +8,14 @@
 *
 **/ 
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 
 function helix(radius, levels, level_dist, vt_dir = "SPI_DOWN", rt_dir = "CT_CLK") = 
     let(
-        is_vt = __is_vector(radius),
-        r1 = is_vt ? radius[0] : radius,
-        r2 = is_vt ? radius[1] : radius,
+        is_flt = __is_float(radius),
+        r1 = is_flt ? radius : radius[0],
+        r2 = is_flt ? radius : radius[1],
         init_r = vt_dir == "SPI_DOWN" ? r2 : r1,
         _frags = __frags(init_r),
         h = level_dist * levels,

src/helix_extrude.scad

@@ -8,7 +8,7 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 
 module helix_extrude(shape_pts, radius, levels, level_dist, 
@@ -22,9 +22,9 @@ module helix_extrude(shape_pts, radius, levels, level_dist,
                 vt[leng - 1 - i]
         ];                         
                          
-    is_vt = __is_vector(radius);
-    r1 = is_vt ? radius[0] : radius;
-    r2 = is_vt ? radius[1] : radius;
+    is_flt = __is_float(radius);
+    r1 = is_flt ? radius : radius[0];
+    r2 = is_flt ? radius : radius[1];
     
     init_r = vt_dir == "SPI_DOWN" ? r2 : r1;
 

src/in_polyline.scad

@@ -0,0 +1,22 @@
+/**
+* in_polyline.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-in_polyline.html
+*
+**/
+
+include <__private__/__to3d.scad>;
+include <__private__/__in_line.scad>;
+
+function _in_polyline_sub(pts, pt, epsilon, iend, i = 0) = 
+    i == iend ? false : (
+        __in_line([pts[i], pts[i + 1]], pt, epsilon) ? true :
+            _in_polyline_sub(pts, pt, epsilon, iend, i + 1)
+    );
+
+function in_polyline(line_pts, pt, epsilon = 0.0001) = 
+    _in_polyline_sub(line_pts, pt, epsilon, len(line_pts) - 1);
+    

src/in_shape.scad

@@ -0,0 +1,56 @@
+/**
+* in_shape.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-in_shape.html
+*
+**/
+
+include <__private__/__to3d.scad>;
+include <__private__/__lines_from.scad>;
+include <__private__/__in_line.scad>;
+
+function _in_shape_in_line_equation(edge, pt) = 
+    let(
+        x1 = edge[0][0],
+        y1 = edge[0][1],
+        x2 = edge[1][0],
+        y2 = edge[1][1],
+        a = (y2 - y1) / (x2 - x1),
+        b = y1 - a * x1
+    )
+    (pt[1] == a * pt[0] + b);
+
+function _in_shape_in_any_edges_sub(edges, leng, pt, i, epsilon) = 
+    leng == i ? false : (
+        __in_line(edges[i], pt, epsilon) ? true : _in_shape_in_any_edges_sub(edges, leng, pt, i + 1, epsilon)
+    );
+
+function _in_shape_in_any_edges(edges, pt, epsilon) = _in_shape_in_any_edges_sub(edges, len(edges), pt, 0, epsilon);
+    
+function _in_shape_interpolate_x(y, p1, p2) = 
+    p1[1] == p2[1] ? p1[0] : (
+        p1[0] + (p2[0] - p1[0]) * (y - p1[1]) / (p2[1] - p1[1])
+    );
+    
+function _in_shape_does_pt_cross(pts, i, j, pt) = 
+    ((pts[i][1] > pt[1]) != (pts[j][1] > pt[1])) &&
+    (pt[0] < _in_shape_interpolate_x(pt[1], pts[i], pts[j]));
+    
+
+function _in_shape_sub(shapt_pts, leng, pt, cond, i, j) =
+    j == leng ? cond : (
+        _in_shape_does_pt_cross(shapt_pts, i, j, pt) ? 
+            _in_shape_sub(shapt_pts, leng, pt, !cond, j, j + 1) :
+            _in_shape_sub(shapt_pts, leng, pt, cond, j, j + 1)
+    );
+ 
+function in_shape(shapt_pts, pt, include_edge = false, epsilon = 0.0001) = 
+    let(
+        leng = len(shapt_pts),
+        edges = __lines_from(points, true)
+    )
+    _in_shape_in_any_edges(edges, pt, epsilon) ? include_edge : 
+    _in_shape_sub(shapt_pts, leng, pt, false, leng - 1, 0);

src/m_rotation.scad

@@ -8,6 +8,9 @@
 *
 **/
 
+include <__private__/__is_float.scad>;
+include <__private__/__to_ang_vect.scad>;
+
 function _q_rotation(a, v) = 
     let(
         half_a = a / 2,
@@ -66,15 +69,8 @@ function _m_zRotation(a) =
         [0, 0, 0, 1]
     ];    
 
-function _to_avect(a) =
-     len(a) == 3 ? a : (
-         len(a) == 2 ? [a[0], a[1], 0] : (
-             len(a) == 1 ? [a[0], 0, 0] : [0, 0, a]
-         ) 
-     );
-
 function _xyz_rotation(a) =
-    let(ang = _to_avect(a))
+    let(ang = __to_ang_vect(a))
     _m_zRotation(ang[2]) * _m_yRotation(ang[1]) * _m_xRotation(ang[0]);
 
 function m_rotation(a, v) = 

src/m_scaling.scad

@@ -8,15 +8,18 @@
 *
 **/
 
-function _to_svect(s) =
-     len(s) == 3 ? s : (
-         len(s) == 2 ? [s[0], s[1], 1] : (
-             len(s) == 1 ? [s[0], 1, 1] : [s, s, s]
-         ) 
+include <__private__/__is_float.scad>;
+
+function _to_3_elems_scaling_vect(s) =
+     let(leng = len(s))
+     leng == 3 ? s : (
+         leng == 2 ? [s[0], s[1], 1] : [s[0], 1, 1]
      );
 
+function _to_scaling_vect(s) = __is_float(s) ? [s, s, s] : _to_3_elems_scaling_vect(s);
+
 function m_scaling(s) = 
-    let(v = _to_svect(s))
+    let(v = _to_scaling_vect(s))
     [
         [v[0], 0, 0, 0],
         [0, v[1], 0, 0],

src/m_translation.scad

@@ -8,15 +8,18 @@
 *
 **/
 
-function _to_tvect(v) =
-     len(v) == 3 ? v : (
-         len(v) == 2 ? [v[0], v[1], 0] : (
-             len(v) == 1 ? [v[0], 0, 0] : [v, 0, 0]
-         ) 
+include <__private__/__is_float.scad>;
+
+function _to_3_elems_translation_vect(v) =
+     let(leng = len(v))
+     leng == 3 ? v : (
+         leng == 2 ? [v[0], v[1], 0] : [v[0], 0, 0]
      );
 
+function _to_translation_vect(v) = __is_float(v) ? [v, 0, 0] : _to_3_elems_translation_vect(v);
+
 function m_translation(v) = 
-    let(vt = _to_tvect(v))
+    let(vt = _to_translation_vect(v))
     [
         [1, 0, 0, vt[0]],
         [0, 1, 0, vt[1]],

src/midpt_smooth.scad

@@ -0,0 +1,21 @@
+/**
+* midpt_smooth.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-midpt_smooth.html
+*
+**/
+
+function _midpt_smooth_sub(points, iend, i, closed = false) = 
+    i == iend ? (
+        closed ? [(points[i] + points[0]) / 2]
+        : []
+    ) : concat([(points[i] + points[i + 1]) / 2], _midpt_smooth_sub(points, iend, i + 1, closed)); 
+
+function midpt_smooth(points, n, closed = false) =
+    let(
+        smoothed = _midpt_smooth_sub(points, len(points) - 1, 0, closed)
+    )
+    n == 1 ? smoothed : midpt_smooth(smoothed, n - 1, closed);

src/path_extrude.scad

@@ -8,7 +8,7 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__to3d.scad>;
 include <__private__/__angy_angz.scad>;
 
@@ -16,150 +16,212 @@ include <__private__/__angy_angz.scad>;
 // For backward compatibility, I directly include m_rotation here.
 include <m_rotation.scad>;
 
-module path_extrude(shape_pts, path_pts, triangles = "SOLID", twist = 0, scale = 1.0, closed = false) {
-    // pre-process parameters
-    function scale_pts(pts, s) = [
-        for(p = pts) [p[0] * s[0], p[1] * s[1], p[2] * s[2]]
-    ];
-    
-    function translate_pts(pts, t) = [
-        for(p = pts) [p[0] + t[0], p[1] + t[1], p[2] + t[2]]
-    ];
-        
-    function rotate_pts(pts, a, v) = [for(p = pts) rotate_p(p, a, v)];
-    
+module path_extrude(shape_pts, path_pts, triangles = "SOLID", twist = 0, scale = 1.0, closed = false, method = "AXIS_ANGLE") {
     sh_pts = len(shape_pts[0]) == 3 ? shape_pts : [for(p = shape_pts) __to3d(p)];
     pth_pts = len(path_pts[0]) == 3 ? path_pts : [for(p = path_pts) __to3d(p)];
         
     len_path_pts = len(pth_pts);
     len_path_pts_minus_one = len_path_pts - 1;
-    twist_step_a = twist / len_path_pts;
-
-    function scale_step() =
-        let(s =  (scale - 1) / len_path_pts_minus_one)
-        [s, s, s];  
-
-    scale_step_vt = __is_vector(scale) ? 
-        [
-            (scale[0] - 1) / len_path_pts_minus_one, 
-            (scale[1] - 1) / len_path_pts_minus_one,
-            scale[2] == undef ? 0 : (scale[2] - 1) / len_path_pts_minus_one
-        ] : scale_step();   
-
-    // get rotation matrice for sections
-
-    function local_ang_vects(j) = 
-        j == 0 ? [] : local_ang_vects_sub(j);
     
-    function local_ang_vects_sub(j) =
-        let(
-            vt0 = pth_pts[j] - pth_pts[j - 1],
-            vt1 = pth_pts[j + 1] - pth_pts[j],
-            a = acos((vt0 * vt1) / (norm(vt0) * norm(vt1))),
-            v = cross(vt0, vt1)
-        )
-        concat([[a, v]], local_ang_vects(j - 1));
 
-    rot_matrice = [
-        for(ang_vect = local_ang_vects(len_path_pts - 2)) 
-            m_rotation(ang_vect[0], ang_vect[1])
-    ];
+    module axis_angle_path_extrude() {
+        twist_step_a = twist / len_path_pts;
 
-    leng_rot_matrice = len(rot_matrice);
-    leng_rot_matrice_minus_one = leng_rot_matrice - 1;
-    leng_rot_matrice_minus_two= leng_rot_matrice - 2;
+        function scale_pts(pts, s) = [
+            for(p = pts) [p[0] * s[0], p[1] * s[1], p[2] * s[2]]
+        ];
+        
+        function translate_pts(pts, t) = [
+            for(p = pts) [p[0] + t[0], p[1] + t[1], p[2] + t[2]]
+        ];
+            
+        function rotate_pts(pts, a, v) = [for(p = pts) rotate_p(p, a, v)];
+        
+        function scale_step() =
+            let(s =  (scale - 1) / len_path_pts_minus_one)
+            [s, s, s];  
 
-    identity_matrix = [
-        [1, 0, 0, 0],
-        [0, 1, 0, 0],
-        [0, 0, 1, 0],
-        [0, 0, 0, 1]
-    ];
+        scale_step_vt = __is_float(scale) ? 
+            scale_step() : 
+            [
+                (scale[0] - 1) / len_path_pts_minus_one, 
+                (scale[1] - 1) / len_path_pts_minus_one,
+                scale[2] == undef ? 0 : (scale[2] - 1) / len_path_pts_minus_one
+            ];   
 
-    function cumulated_rot_matrice(i) = 
-        leng_rot_matrice == 0 ? [identity_matrix] : (
-             leng_rot_matrice == 1 ? [rot_matrice[0], identity_matrix] : 
-                (
-                    i == leng_rot_matrice_minus_two ? 
-                    [
-                        rot_matrice[leng_rot_matrice_minus_one], 
-                        rot_matrice[leng_rot_matrice_minus_two] * rot_matrice[leng_rot_matrice_minus_one]
-                    ] 
-                    : cumulated_rot_matrice_sub(i))
-        );
+        // get rotation matrice for sections
 
-    function cumulated_rot_matrice_sub(i) = 
-        let(
-            matrice = cumulated_rot_matrice(i + 1),
-            curr_matrix = rot_matrice[i],
-            prev_matrix = matrice[len(matrice) - 1]
-        )
-        concat(matrice, [curr_matrix * prev_matrix]);
+        function local_ang_vects(j) = 
+            j == 0 ? [] : local_ang_vects_sub(j);
+        
+        function local_ang_vects_sub(j) =
+            let(
+                vt0 = pth_pts[j] - pth_pts[j - 1],
+                vt1 = pth_pts[j + 1] - pth_pts[j],
+                a = acos((vt0 * vt1) / (norm(vt0) * norm(vt1))),
+                v = cross(vt0, vt1)
+            )
+            concat([[a, v]], local_ang_vects(j - 1));
 
-    cumu_rot_matrice = cumulated_rot_matrice(0);
+        rot_matrice = [
+            for(ang_vect = local_ang_vects(len_path_pts - 2)) 
+                m_rotation(ang_vect[0], ang_vect[1])
+        ];
 
-    // get all sections
+        leng_rot_matrice = len(rot_matrice);
+        leng_rot_matrice_minus_one = leng_rot_matrice - 1;
+        leng_rot_matrice_minus_two= leng_rot_matrice - 2;
 
-    function init_section(a, s) =
-        let(angleyz = __angy_angz(pth_pts[0], pth_pts[1]))
-        rotate_pts(
+        identity_matrix = [
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, 1, 0],
+            [0, 0, 0, 1]
+        ];
+
+        function cumulated_rot_matrice(i) = 
+            leng_rot_matrice == 0 ? [identity_matrix] : (
+                leng_rot_matrice == 1 ? [rot_matrice[0], identity_matrix] : 
+                    (
+                        i == leng_rot_matrice_minus_two ? 
+                        [
+                            rot_matrice[leng_rot_matrice_minus_one], 
+                            rot_matrice[leng_rot_matrice_minus_two] * rot_matrice[leng_rot_matrice_minus_one]
+                        ] 
+                        : cumulated_rot_matrice_sub(i))
+            );
+
+        function cumulated_rot_matrice_sub(i) = 
+            let(
+                matrice = cumulated_rot_matrice(i + 1),
+                curr_matrix = rot_matrice[i],
+                prev_matrix = matrice[len(matrice) - 1]
+            )
+            concat(matrice, [curr_matrix * prev_matrix]);
+
+        cumu_rot_matrice = cumulated_rot_matrice(0);
+
+        // get all sections
+
+        function init_section(a, s) =
+            let(angleyz = __angy_angz(pth_pts[0], pth_pts[1]))
             rotate_pts(
                 rotate_pts(
-                    scale_pts(sh_pts, s), a
-                ), [90, 0, -90]
-            ), [0, -angleyz[0], angleyz[1]]
-        );
+                    rotate_pts(
+                        scale_pts(sh_pts, s), a
+                    ), [90, 0, -90]
+                ), [0, -angleyz[0], angleyz[1]]
+            );
+            
+        function local_rotate_section(j, init_a, init_s) =
+            j == 0 ? 
+                init_section(init_a, init_s) : 
+                local_rotate_section_sub(j, init_a, init_s);
         
-    function local_rotate_section(j, init_a, init_s) =
-        j == 0 ? 
-            init_section(init_a, init_s) : 
-            local_rotate_section_sub(j, init_a, init_s);
-    
-    function local_rotate_section_sub(j, init_a, init_s) = 
-        let(
-            vt0 = pth_pts[j] - pth_pts[j - 1],
-            vt1 = pth_pts[j + 1] - pth_pts[j],
-            ms = cumu_rot_matrice[j - 1]
-        )
-        [
-            for(p = init_section(init_a, init_s)) 
-                [
-                    [ms[0][0], ms[0][1], ms[0][2]] * p,
-                    [ms[1][0], ms[1][1], ms[1][2]] * p,
-                    [ms[2][0], ms[2][1], ms[2][2]] * p
-                ]
-        ];        
+        function local_rotate_section_sub(j, init_a, init_s) = 
+            let(
+                vt0 = pth_pts[j] - pth_pts[j - 1],
+                vt1 = pth_pts[j + 1] - pth_pts[j],
+                ms = cumu_rot_matrice[j - 1]
+            )
+            [
+                for(p = init_section(init_a, init_s)) 
+                    [
+                        [ms[0][0], ms[0][1], ms[0][2]] * p,
+                        [ms[1][0], ms[1][1], ms[1][2]] * p,
+                        [ms[2][0], ms[2][1], ms[2][2]] * p
+                    ]
+            ];        
 
-    function sections() =
-        let(
-            fst_section = 
-                translate_pts(local_rotate_section(0, 0, [1, 1, 1]), pth_pts[0]),
-            remain_sections = [
-                for(i = [0:len_path_pts - 2]) 
-                
-                    translate_pts(
-                        local_rotate_section(i, i * twist_step_a, [1, 1, 1] + scale_step_vt * i),
-                        pth_pts[i + 1]
-                    )
-                
+        function sections() =
+            let(
+                fst_section = 
+                    translate_pts(local_rotate_section(0, 0, [1, 1, 1]), pth_pts[0]),
+                remain_sections = [
+                    for(i = [0:len_path_pts - 2]) 
                     
-            ]
-        ) concat([fst_section], remain_sections);
+                        translate_pts(
+                            local_rotate_section(i, i * twist_step_a, [1, 1, 1] + scale_step_vt * i),
+                            pth_pts[i + 1]
+                        )
+                    
+                        
+                ]
+            ) concat([fst_section], remain_sections);
+        
+        sects = sections();
+
+        function calculated_sections() =
+            closed && pth_pts[0] == pth_pts[len_path_pts_minus_one] ?
+                concat(sects, [sects[0]]) : // round-robin
+                sects;
+        
+        polysections(
+            calculated_sections(),
+            triangles = triangles
+        );   
+
+        // hook for testing
+        test_path_extrude(sects);        
+    }
+
+    module euler_angle_path_extrude() {
+        scale_step_vt = __is_float(scale) ? 
+            [(scale - 1) / len_path_pts_minus_one, (scale - 1) / len_path_pts_minus_one] : 
+            [(scale[0] - 1) / len_path_pts_minus_one, (scale[1] - 1) / len_path_pts_minus_one];
+
+        scale_step_x = scale_step_vt[0];
+        scale_step_y = scale_step_vt[1];
+        twist_step = twist / len_path_pts_minus_one;
+
+        function section(p1, p2, i) = 
+            let(
+                length = norm(p1 - p2),
+                angy_angz = __angy_angz(p1, p2),
+                ay = -angy_angz[0],
+                az = angy_angz[1]
+            )
+            [
+                for(p = sh_pts) 
+                    let(scaled_p = [p[0] * (1 + scale_step_x * i), p[1] * (1 + scale_step_y * i), p[2]])
+                    rotate_p(
+                        rotate_p(
+                            rotate_p(scaled_p, twist_step * i), [90, 0, -90]
+                        ) + [i == 0 ? 0 : length, 0, 0], 
+                        [0, ay, az]
+                    ) + p1
+            ];
+        
+        function path_extrude_inner(index) =
+        index == len_path_pts ? [] :
+            concat(
+                [section(pth_pts[index - 1], pth_pts[index],  index)],
+                path_extrude_inner(index + 1)
+            );
+
+        function calculated_sections() =
+            let(sections = path_extrude_inner(1))
+            closed && pth_pts[0] == pth_pts[len_path_pts_minus_one] ?
+                concat(sections, [sections[0]]) : // round-robin
+                concat([section(pth_pts[0], pth_pts[1], 0)], sections);
     
-    sects = sections();
+    sections = calculated_sections();
 
-    function calculated_sections() =
-        closed && pth_pts[0] == pth_pts[len_path_pts_minus_one] ?
-            concat(sects, [sects[0]]) : // round-robin
-            sects;
-     
-     polysections(
-        calculated_sections(),
-        triangles = triangles
-    );   
+        polysections(
+            sections,
+            triangles = triangles
+        );   
 
-    // hook for testing
-    test_path_extrude(sects);
+        // hook for testing
+        test_path_extrude(sections);
+    }
+
+    if(method == "AXIS_ANGLE") {
+        axis_angle_path_extrude();
+    }
+    else if(method == "EULER_ANGLE") {
+        euler_angle_path_extrude();
+    } 
 }
 
 

src/path_scaling_sections.scad

@@ -0,0 +1,30 @@
+/**
+* path_scaling_sections.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-path_scaling_sections.html
+*
+**/
+
+include <__private__/__reverse.scad>;
+
+function path_scaling_sections(shape_pts, edge_path) = 
+    let(
+        start_point = edge_path[0],
+        base_leng = norm(start_point),
+        scaling_matrice = [
+            for(p = edge_path) 
+            let(s = norm([p[0], p[1], 0]) / base_leng)
+            m_scaling([s, s, 1])
+        ]
+    )
+    __reverse([
+        for(i = [0:len(edge_path) - 1])
+        [
+            for(p = shape_pts) 
+            let(scaled_p = scaling_matrice[i] * [p[0], p[1], edge_path[i][2], 1])
+            [scaled_p[0], scaled_p[1], scaled_p[2]]
+        ]
+    ]);

src/pie.scad

@@ -9,7 +9,7 @@
 **/
 
 include <__private__/__frags.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__ra_to_xy.scad>;
 include <__private__/__shape_pie.scad>;
  

src/ring_extrude.scad

@@ -19,7 +19,7 @@ module ring_extrude(shape_pts, radius, angle = 360, twist = 0, scale = 1.0, tria
     } else {
         a_step = 360 / __frags(radius);
 
-        angles = __is_vector(angle) ? angle : [0, angle];
+        angles = __is_float(angle) ? [0, angle] : angle;
 
         m = floor(angles[0] / a_step) + 1;
         n = floor(angles[1] / a_step);

src/rotate_p.scad

@@ -10,7 +10,8 @@
 
 include <__private__/__to2d.scad>;
 include <__private__/__to3d.scad>;
-
+include <__private__/__is_float.scad>;
+include <__private__/__to_ang_vect.scad>;
 
 function _q_rotate_p_3d(p, a, v) = 
     let(
@@ -69,15 +70,8 @@ function _rotz(pt, a) =
 function _rotate_p_3d(point, a) =
     _rotz(_roty(_rotx(point, a[0]), a[1]), a[2]);
 
-function _to_avect(a) =
-     len(a) == 3 ? a : (
-         len(a) == 2 ? [a[0], a[1], 0] : (
-             len(a) == 1 ? [a[0], 0, 0] : [0, 0, a]
-         ) 
-     );
-
 function _rotate_p(p, a) =
-    let(angle = _to_avect(a))
+    let(angle = __to_ang_vect(a))
     len(p) == 3 ? 
         _rotate_p_3d(p, angle) :
         __to2d(

src/rounded_cube.scad

@@ -8,15 +8,15 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 include <__private__/__nearest_multiple_of_4.scad>;
 
 module rounded_cube(size, corner_r, center = false) {
-    is_vt = __is_vector(size);
-    x = is_vt ? size[0] : size;
-    y = is_vt ? size[1] : size;
-    z = is_vt ? size[2] : size;
+    is_flt = __is_float(size);
+    x = is_flt ? size : size[0];
+    y = is_flt ? size : size[1];
+    z = is_flt ? size : size[2];
 
     corner_frags = __nearest_multiple_of_4(__frags(corner_r));
     edge_d = corner_r * cos(180 / corner_frags);

src/rounded_cylinder.scad

@@ -8,7 +8,7 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 include <__private__/__pie_for_rounding.scad>;
 include <__private__/__half_trapezium.scad>;

src/rounded_extrude.scad

@@ -9,13 +9,13 @@
 **/
 
 include <__private__/__frags.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 
 module rounded_extrude(size, round_r, angle = 90, twist = 0, convexity = 10) {
 
-    is_vt = __is_vector(size);
-    x = is_vt ? size[0] : size;
-    y = is_vt ? size[1] : size;
+    is_flt = __is_float(size);
+    x = is_flt ? size : size[0];
+    y = is_flt ? size : size[1];
     
     q_corner_frags = __frags(round_r) / 4;
     

src/rounded_square.scad

@@ -8,16 +8,16 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 include <__private__/__pie_for_rounding.scad>;
 include <__private__/__half_trapezium.scad>;
 include <__private__/__trapezium.scad>;
 
 module rounded_square(size, corner_r, center = false) {
-    is_vt = __is_vector(size);
-    x = is_vt ? size[0] : size;
-    y = is_vt ? size[1] : size;       
+    is_flt = __is_float(size);
+    x = is_flt ? size : size[0];
+    y = is_flt ? size : size[1];       
     
     position = center ? [0, 0] : [x / 2, y / 2];
     points = __trapezium(

src/shape_arc.scad

@@ -9,9 +9,10 @@
 **/ 
 
 include <__private__/__frags.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__ra_to_xy.scad>;
 include <__private__/__shape_arc.scad>;
+include <__private__/__edge_r.scad>
 
 function shape_arc(radius, angle, width, width_mode = "LINE_CROSS") =
     __shape_arc(radius, angle, width, width_mode);

src/shape_pie.scad

@@ -9,7 +9,7 @@
 **/
 
 include <__private__/__frags.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__ra_to_xy.scad>;
 include <__private__/__shape_pie.scad>;
 

src/shape_square.scad

@@ -8,17 +8,17 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 include <__private__/__pie_for_rounding.scad>;
 include <__private__/__half_trapezium.scad>;
 include <__private__/__trapezium.scad>;
-
+ 
 function shape_square(size, corner_r = 0) = 
     let(
-        is_vt = __is_vector(size),
-        x = is_vt ? size[0] : size,
-        y = is_vt ? size[1] : size        
+        is_flt = __is_float(size),
+        x = is_flt ? size : size[0],
+        y = is_flt ? size : size[1]        
     )
     __trapezium(
         length = x, 

src/shape_trapezium.scad

@@ -8,7 +8,7 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__frags.scad>;
 include <__private__/__pie_for_rounding.scad>;
 include <__private__/__half_trapezium.scad>;

src/starburst.scad

@@ -0,0 +1,42 @@
+/**
+* starburst.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-starburst.html
+*
+**/
+
+module starburst(r1, r2, n, height) {
+    a = 180 / n;
+
+    p0 = [0, 0, 0];
+    p1 = [r2 * cos(a), r2 * sin(a), 0];
+    p2 = [r1, 0, 0];
+    p3 = [0, 0, height];
+
+    module half_burst() {
+        polyhedron(points = [p0, p1, p2, p3], 
+            faces = [
+                [0, 2, 1],
+                [0, 1, 3],
+                [0, 3, 2], 
+                [2, 1, 3]
+            ]
+        );
+    }
+
+    module burst() {
+        hull() {
+            half_burst();
+            mirror([0, 1,0]) half_burst();
+        }
+    }
+
+    union() {
+        for(i = [0 : n - 1]) {
+            rotate(2 * a * i) burst();
+        }
+    }
+}

src/torus_knot.scad

@@ -0,0 +1,21 @@
+/**
+* torus_knot.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-torus_knot.html
+*
+**/
+
+function torus_knot(p, q, phi_step) = [
+    for(phi = [0:phi_step:6.28318])
+    let(
+        degree = phi * 180 / 3.14159,
+        r = cos(q * degree) + 2,
+        x = r * cos(p * degree),
+        y = r * sin(p * degree),
+        z = -sin(q * degree)
+    )
+    [x, y, z]
+];

src/triangulate.scad

@@ -0,0 +1,103 @@
+/**
+* triangulate.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-triangulate.html
+*
+**/
+
+function _triangulate_in_triangle(p0, p1, p2, p) =
+    let(
+        v0 = p0 - p,
+        v1 = p1 - p,
+        v2 = p2 - p,
+        c0 = cross(v0, v1),
+        c1 = cross(v1, v2),
+        c2 = cross(v2, v0)
+    )
+    (c0 > 0 && c1 > 0 && c2 > 0) || (c0 < 0 && c1 < 0 && c2 < 0);
+
+function _triangulate_snipable(shape_pts, u, v, w, n, indices, epsilon = 0.0001) =
+    let(
+        a = shape_pts[indices[u]],
+        b = shape_pts[indices[v]],
+        c = shape_pts[indices[w]],
+        ax = a[0],
+        ay = a[1],
+        bx = b[0],
+        by = b[1],
+        cx = c[0],
+        cy = c[1]
+    )
+    epsilon > (((bx - ax) * (cy - ay)) - ((by - ay) * (cx - ax))) ? 
+        false : _triangulate_snipable_sub(shape_pts, n, u, v, w, a, b, c, indices);
+    
+function _triangulate_snipable_sub(shape_pts, n, u, v, w, a, b, c, indices, p = 0) = 
+    p == n ? true : (
+        ((p == u) || (p == v) || (p == w)) ? _triangulate_snipable_sub(shape_pts, n, u, v, w, a, b, c, indices, p + 1) : (
+            _triangulate_in_triangle(a, b, c, shape_pts[indices[p]]) ? 
+                false : _triangulate_snipable_sub(shape_pts, n, u, v, w, a, b, c, indices, p + 1)
+        )
+    );
+
+// remove the elem at idx v from indices  
+function _triangulate_remove_v(indices, v, num_of_vertices) = 
+    let(
+        nv_minuns_one = num_of_vertices - 1
+    )
+    v == 0 ? [for(i = [1:nv_minuns_one]) indices[i]] : (
+        v == nv_minuns_one ? [for(i = [0:v - 1]) indices[i]] : concat(
+            [for(i = [0:v - 1]) indices[i]], 
+            [for(i = [v + 1:nv_minuns_one]) indices[i]]
+        )
+    );
+    
+function _triangulate_zero_or_value(num_of_vertices, value) = 
+    num_of_vertices <= value ? 0 : value;
+
+
+function _triangulate_real_triangulate_sub(shape_pts, collector, indices, v, num_of_vertices, count, epsilon) = 
+    let(
+        // idxes of three consecutive vertices
+        u = _triangulate_zero_or_value(num_of_vertices, v),     
+        vi = _triangulate_zero_or_value(num_of_vertices, u + 1), 
+        w = _triangulate_zero_or_value(num_of_vertices, vi + 1) 
+    )
+    _triangulate_snipable(shape_pts, u, vi, w, num_of_vertices, indices, epsilon) ? 
+        _triangulate_snip(shape_pts, collector, indices, u, vi, w, num_of_vertices, count, epsilon) : 
+        _triangulate_real_triangulate(shape_pts, collector, indices, vi, num_of_vertices, count, epsilon);
+        
+function _triangulate_snip(shape_pts, collector, indices, u, v, w, num_of_vertices, count, epsilon) = 
+    let(
+        a = indices[u],
+        b = indices[v],
+        c = indices[w],
+        new_nv = num_of_vertices - 1
+    )
+    _triangulate_real_triangulate( 
+        shape_pts, 
+        concat(collector, [[a, b, c]]),  
+        _triangulate_remove_v(indices, v, num_of_vertices),
+        v, 
+        new_nv,
+        2 * new_nv,
+        epsilon
+    );
+
+function _triangulate_real_triangulate(shape_pts, collector, indices, v, num_of_vertices, count, epsilon) = 
+    count <= 0 ? [] : (
+        num_of_vertices == 2 ? 
+            collector : _triangulate_real_triangulate_sub(shape_pts, collector, indices, v, num_of_vertices, count - 1,  epsilon)
+    );
+    
+function triangulate(shape_pts,  epsilon = 0.0001) =
+    let(
+        num_of_vertices = len(shape_pts),
+        v = num_of_vertices - 1,
+        indices = [for(vi = [0:v]) vi],
+        count = 2 * num_of_vertices        
+    )
+    num_of_vertices < 3 ? [] : _triangulate_real_triangulate(shape_pts, [], indices, v, num_of_vertices, count, epsilon);
+       

src/trim_shape.scad

@@ -0,0 +1,59 @@
+/**
+* trim_shape.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-trim_shape.html
+*
+**/
+
+include <__private__/__to3d.scad>;
+include <__private__/__line_intersection.scad>;
+include <__private__/__in_line.scad>;
+include <__private__/__lines_from.scad>;
+
+function _trim_shape_any_intersection_sub(lines, line, lines_leng, i, epsilon) =
+    let(
+        p = __line_intersection(lines[i], line, epsilon)
+    )
+    (p != [] && __in_line(line, p, epsilon) && __in_line(lines[i], p, epsilon)) ? [i, p] : _trim_shape_any_intersection(lines, line, lines_leng, i + 1, epsilon);
+
+// return [idx, [x, y]] or []
+function _trim_shape_any_intersection(lines, line, lines_leng, i, epsilon) =
+    i == lines_leng ? [] : _trim_shape_any_intersection_sub(lines, line, lines_leng, i, epsilon);
+
+function _trim_sub(lines, leng, epsilon) = 
+    let(
+        current_line = lines[0],
+        next_line = lines[1],
+        lines_from_next = [for(j = [1 : leng - 1]) lines[j]],
+        lines_from_next2 = [for(j = [2 : leng - 1]) lines[j]],
+        current_p = current_line[0],
+        leng_lines_from_next2 = len(lines_from_next2),
+        inter_p = _trim_shape_any_intersection(lines_from_next2, current_line, leng_lines_from_next2, 0, epsilon)
+    )
+    // no intersecting pt, collect current_p and trim remain lines
+    inter_p == [] ? (concat([current_p], _trim_shape_trim_lines(lines_from_next, epsilon))) : (
+        // collect current_p, intersecting pt and the last pt
+        (leng == 3 || (inter_p[0] == (leng_lines_from_next2 - 1))) ? [current_p, inter_p[1], lines[leng - 1]] : (
+            // collect current_p, intersecting pt and trim remain lines
+            concat([current_p, inter_p[1]], 
+                _trim_shape_trim_lines([for(i = [inter_p[0] + 1 : leng_lines_from_next2 - 1]) lines_from_next2[i]], epsilon)
+            )
+        )
+    );
+    
+function _trim_shape_trim_lines(lines, epsilon) = 
+    let(leng = len(lines))
+    leng > 2 ? _trim_sub(lines, leng, epsilon) : _trim_shape_collect_pts_from(lines, leng);
+
+function _trim_shape_collect_pts_from(lines, leng) = 
+    concat([for(line = lines) line[0]], [lines[leng - 1][1]]);
+
+function trim_shape(shape_pts, from, to, epsilon = 0.0001) = 
+    let(
+        pts = [for(i = [from:to]) shape_pts[i]],
+        trimmed = _trim_shape_trim_lines(__lines_from(pts), epsilon)
+    )
+    len(shape_pts) == len(trimmed) ? trimmed : trim_shape(trimmed, 0, len(trimmed) - 1, epsilon);

src/voronoi2d.scad

@@ -0,0 +1,46 @@
+/**
+* voronoi2d.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-voronoi2d.html
+*
+**/
+
+module voronoi2d(points, spacing = 1, r = 0, delta = 0, chamfer = false, region_type = "square") {
+    xs = [for(p = points) p[0]];
+    ys = [for(p = points) abs(p[1])];
+
+    region_size = max([(max(xs) -  min(xs) / 2), (max(ys) -  min(ys)) / 2]);    
+    half_region_size = 0.5 * region_size; 
+    offset_leng = spacing * 0.5 + half_region_size;
+
+    function normalize(v) = v / norm(v);
+    
+    module region(pt) {
+        intersection_for(p = points) {
+            if(pt != p) {
+                v = p - pt;
+                translate((pt + p) / 2 - normalize(v) * offset_leng)
+                    rotate(atan2(v[1], v[0])) 
+                    if(region_type == "square") {
+                        square(region_size, center = true);
+                    }
+                    else if(region_type == "circle") {
+                        circle(region_size / 2);
+                    }
+                    
+            }
+        }
+    }    
+    
+    for(p = points) {	
+        if(r != 0) {        
+            offset(r) region(p);
+        }
+        else {
+            offset(delta = delta, chamfer = chamfer) region(p);
+        }
+    }
+}

src/voronoi3d.scad

@@ -0,0 +1,42 @@
+/**
+* voronoi3d.scad
+*
+* @copyright Justin Lin, 2019
+* @license https://opensource.org/licenses/lgpl-3.0.html
+*
+* @see https://openhome.cc/eGossip/OpenSCAD/lib-voronoi3d.html
+*
+**/
+
+include <__private__/__angy_angz.scad>;
+
+// slow but workable
+
+module voronoi3d(points, space_size = "auto", spacing = 1) {
+    xs = [for(p = points) p[0]];
+    ys = [for(p = points) abs(p[1])];
+    zs = [for(p = points) abs(p[2])];
+
+    space_size = max([max(xs) -  min(xs), max(ys) -  min(ys), max(zs) -  min(zs)]);    
+    half_space_size = 0.5 * space_size; 
+    offset_leng = spacing * 0.5 + half_space_size;
+
+    function normalize(v) = v / norm(v);
+    
+    module space(pt) {
+        intersection_for(p = points) {
+            if(pt != p) {
+                v = p - pt;
+                ryz = __angy_angz(p, pt);
+
+                translate((pt + p) / 2 - normalize(v) * offset_leng)
+                    rotate([0, -ryz[0], ryz[1]]) 
+                    cube([space_size, space_size * 2, space_size * 2], center = true); 
+            }
+        }
+    }    
+    
+    for(p = points) {	
+        space(p);
+    }
+}

test/test_hull_polyline2d.scad

@@ -14,8 +14,6 @@ module test_hull_polyline2d() {
         assertEqualPoint(points[index], point2);
         assertEqual(line_width, radius * 2);    
     } 
-    
-    $fn = 4;
 
     hull_polyline2d(
         points = [[1, 2], [-5, -4], [-5, 3], [5, 5]],