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)
 
@@ -39,6 +39,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 +50,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 +66,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 +84,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 +107,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)
@@ -104,7 +116,6 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [sphere_spiral_extrude](https://openhome.cc/eGossip/OpenSCAD/lib-sphere_spiral_extrude.html)
 
 - Matrix
-	- [m_multiply](https://openhome.cc/eGossip/OpenSCAD/lib-m_multiply.html)
 	- [m_cumulate](https://openhome.cc/eGossip/OpenSCAD/lib-m_cumulate.html)	
 	- [m_translation](https://openhome.cc/eGossip/OpenSCAD/lib-m_translation.html)
 	- [m_rotation](https://openhome.cc/eGossip/OpenSCAD/lib-m_rotation.html)

RELEASE.md

@@ -1,3 +1,43 @@
+> Version numbers are based on [Semantic Versioning](https://semver.org/).
+
+# 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.
+  - The `path_pts` parameter of `path_extrude` accepts two or three points. 
+  - The `points` parameter of `along_with` accepts two or three points. 
+
+- Others
+  -  OpenSCAD has built-in matrix multiplication so `m_multiply` is not necessary.
+
 # v1.1
 - New matrix functions:
   - [m_multiply](https://openhome.cc/eGossip/OpenSCAD/lib-m_multiply.html)
@@ -14,7 +54,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)  

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
 

docs/lib-arc_path.md

@@ -15,8 +15,8 @@ Creates an arc path. You can pass a 2 element vector to define the central angle
     include <hull_polyline2d.scad>;
     
     $fn = 24;
-    points = arc_path(radius = 20, angle = [45, 290], width = 2);
-    hull_polyline2d(points);
+    points = arc_path(radius = 20, angle = [45, 290]);
+    hull_polyline2d(points, width = 2);
 
 ![arc_path](images/lib-arc_path-1.JPG)
 
@@ -24,8 +24,8 @@ Creates an arc path. You can pass a 2 element vector to define the central angle
     include <hull_polyline2d.scad>;
     
     $fn = 24;
-    points = arc_path(radius = 20, angle = 135, width = 2);
-    hull_polyline2d(points);
+    points = arc_path(radius = 20, angle = 135);
+    hull_polyline2d(points, width = 2);
 
 ![arc_path](images/lib-arc_path-2.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

@@ -1,23 +0,0 @@
-# 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_cyclicpolygon.md

@@ -6,6 +6,7 @@ Returns shape points of a regular cyclic polygon. They can be used with xxx_extr
 
 - `sides` : The radius of the circle.
 - `circle_r` : The radius of the circumcircle.
+- `corner_r` : The radius of the circle at a corner.
 - `$fa`, `$fs`, `$fn` : Check [the circle module](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Using_the_2D_Subsystem#circle) for more details.
 
 ## Examples

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__/__m_multiply.scad

@@ -1,13 +0,0 @@
-function __m_multiply(ma, mb) = 
-    let(
-        c1 = [mb[0][0], mb[1][0], mb[2][0], mb[3][0]],
-        c2 = [mb[0][1], mb[1][1], mb[2][1], mb[3][1]],
-        c3 = [mb[0][2], mb[1][2], mb[2][2], mb[3][2]],
-        c4 = [mb[0][3], mb[1][3], mb[2][3], mb[3][3]]
-    )
-    [
-        [ma[0] * c1, ma[0] * c2, ma[0] * c3, ma[0] * c4],
-        [ma[1] * c1, ma[1] * c2, ma[1] * c3, ma[1] * c4],
-        [ma[2] * c1, ma[2] * c2, ma[2] * c3, ma[2] * c4],
-        [ma[3] * c1, ma[3] * c2, ma[3] * c3, ma[3] * c4]
-    ];

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

@@ -7,70 +7,91 @@
 * @see https://openhome.cc/eGossip/OpenSCAD/lib-along_with.html
 *
 **/ 
-
+ 
 include <__private__/__angy_angz.scad>;
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__to3d.scad>;
-include <__private__/__m_multiply.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],
+        [0, 0, 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,
             leng_rot_matrice_minus_two = leng_rot_matrice - 2
         )
-        i == leng_rot_matrice - 2 ? 
+        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], 
-                   __m_multiply(rot_matrice[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, 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]
         )
-        concat(matrice, [__m_multiply(curr_matrix, prev_matrix)]);
+        concat(matrice, [curr_matrix * prev_matrix]);
 
     // 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) 
@@ -78,7 +99,7 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
                             children(0);
     }
 
-    module align_with_pts_init(a, s) {
+    module axis_angle_align_with_pts_init(a, s) {
         angleyz = __angy_angz(points[0], points[1]);
         rotate([0, -angleyz[0], angleyz[1]])
             rotate([90, 0, -90])
@@ -87,51 +108,104 @@ 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) {
-        if($children == 1) { 
-            for(i = [0:leng_points_minus_one]) {
-                align_with_pts_angles(i) children(0);
+    // <<< 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") {
+        if($children == 1) { 
+            angs = angles_defined ? angles : euler_angle_path_angles($children);
+            
+            for(i = [0:leng_points_minus_one]) {
+                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(pts) - 1);
+    

src/in_shape.scad

@@ -0,0 +1,55 @@
+/**
+* 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__/__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_cumulate.scad

@@ -8,12 +8,10 @@
 *
 **/
 
-include <__private__/__m_multiply.scad>;
-
 function _m_cumulate(matrice, i) = 
     i == len(matrice) - 2 ?
-        __m_multiply(matrice[i], matrice[i + 1]) :
-        __m_multiply(matrice[i], _m_cumulate(matrice, i + 1));
+        matrice[i] * matrice[i + 1] :
+        matrice[i] * _m_cumulate(matrice, i + 1);
 
 function m_cumulate(matrice) = 
     len(matrice) == 1 ? matrice[0] : _m_cumulate(matrice, 0);

src/m_multiply.scad

@@ -1,13 +0,0 @@
-/**
-* m_multiply.scad
-*
-* @copyright Justin Lin, 2019
-* @license https://opensource.org/licenses/lgpl-3.0.html
-*
-* @see https://openhome.cc/eGossip/OpenSCAD/lib-m_multiply.html
-*
-**/
-
-include <__private__/__m_multiply.scad>;
-
-function m_multiply(ma, mb) = __m_multiply(ma, mb);

src/m_rotation.scad

@@ -8,7 +8,8 @@
 *
 **/
 
-include <__private__/__m_multiply.scad>;
+include <__private__/__is_float.scad>;
+include <__private__/__to_ang_vect.scad>;
 
 function _q_rotation(a, v) = 
     let(
@@ -68,20 +69,14 @@ 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))
-    __m_multiply(
-        _m_zRotation(ang[2]), __m_multiply(
-            _m_yRotation(ang[1]), _m_xRotation(ang[0])
-        )
-    );
+    let(ang = __to_ang_vect(a))
+    _m_zRotation(ang[2]) * _m_yRotation(ang[1]) * _m_xRotation(ang[0]);
 
 function m_rotation(a, v) = 
-    v == undef ? _xyz_rotation(a) : _q_rotation(a, v);
+    (a == 0 || a == [0, 0, 0] || a == [0] || a == [0, 0]) ? [
+        [1, 0, 0, 0],
+        [0, 1, 0, 0],
+        [0, 0, 1, 0],
+        [0, 0, 0, 1]
+    ] : (v == undef ? _xyz_rotation(a) : _q_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_shearing.scad

@@ -8,7 +8,6 @@
 *
 **/
 
-include <__private__/__m_multiply.scad>;
 include <__private__/__m_shearing.scad>;
 
 function m_shearing(sx = [0, 0], sy = [0, 0], sz = [0, 0]) = __m_shearing(sx, sy, sz);

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,148 +8,220 @@
 *
 **/
 
-include <__private__/__is_vector.scad>;
+include <__private__/__is_float.scad>;
 include <__private__/__to3d.scad>;
 include <__private__/__angy_angz.scad>;
-include <__private__/__m_multiply.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 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];  
 
-    function cumulated_rot_matrice(i) = 
-        i == leng_rot_matrice - 2 ? 
-               [
-                   rot_matrice[leng_rot_matrice_minus_one], 
-                   __m_multiply(rot_matrice[leng_rot_matrice_minus_two], rot_matrice[leng_rot_matrice_minus_one])
-               ] 
-               : cumulated_rot_matrice_sub(i);
+        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_sub(i) = 
-        let(
-            matrice = cumulated_rot_matrice(i + 1),
-            curr_matrix = rot_matrice[i],
-            prev_matrix = matrice[len(matrice) - 1]
-        )
-        concat(matrice, [__m_multiply(curr_matrix, prev_matrix)]);
+        // get rotation matrice for sections
 
-    cumu_rot_matrice = cumulated_rot_matrice(0);
+        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));
 
-    // get all sections
+        rot_matrice = [
+            for(ang_vect = local_ang_vects(len_path_pts - 2)) 
+                m_rotation(ang_vect[0], ang_vect[1])
+        ];
 
-    function init_section(a, s) =
-        let(angleyz = __angy_angz(pth_pts[0], pth_pts[1]))
-        rotate_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;
+
+        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/shear.scad

@@ -8,7 +8,6 @@
 *
 **/
 
-include <__private__/__m_multiply.scad>;
 include <__private__/__m_shearing.scad>;
 
 module shear(sx = [0, 0], sy = [0, 0], sz = [0, 0]) {

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);