update RELEASE

README.md

@@ -1,6 +1,6 @@
-# dotSCAD 1.2
+# dotSCAD 1.3
 
-> Helpful modules and functions when playing OpenSCAD. Based on OpenSCAD 2015.03.
+> Reduce the burden of 3D modeling in mathematics. Compatible with OpenSCAD 2015.03 or laters.
 
 ![dotSCAD](WhirlingTaiwan.JPG)
 
@@ -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)
@@ -50,6 +53,7 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [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)
@@ -66,6 +70,11 @@ Too many dependencies? Because OpenSCAD doesn't provide namespace management, I
 	- [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)
@@ -100,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,9 +1,41 @@
 > Version numbers are based on [Semantic Versioning](https://semver.org/).
 
-# v1.2.1
+# v1.3.3
 - Bugfixes
-   - Fixed CCW faces when using `path_scaling_sections`.
+  - `in_shape`: Wrong variable name.
 
+# 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)
@@ -39,7 +71,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-arc.md

@@ -6,6 +6,7 @@ Creates an arc. You can pass a 2 element vector to define the central angle. Its
 
 - `radius` : The radius of the circle.
 - `angle` : A single value or a 2 element vector which defines the central angle. The first element of the vector is the beginning angle in degrees, and the second element is the ending angle.
+- `width` : The width of the arc.
 - `width_mode` : The default value is `"LINE_CROSS"`. The arc line will move outward by `width / 2` and inward by `width / 2`. If it's `"LINE_OUTWARD"`, The arc line moves outward by `width`. The `"LINE_INWARD"` moves the arc line inward by `width`.
 - `$fa`, `$fs`, `$fn` : Check [the circle module](https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Using_the_2D_Subsystem#circle) for more details.
 

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

@@ -8,6 +8,7 @@ Move 2D outlines outward or inward by a given amount. Each point of the offsette
 
 - `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
 

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-m_multiply.md

@@ -1,25 +0,0 @@
-# m_multiply
-
-Multiply two 4x4 transformation matrice.
-
-**Since:** 1.1
-
-## Parameters
-
-- `ma`, `mb` : Two 4x4 transformation matrice.
-
-## Examples
-
-	include <m_multiply.scad>;
-	include <m_scaling.scad>;
-	include <m_rotation.scad>;
-
-	ma = m_scaling([0.5, 1, 2]);
-	mb = m_rotation([0, 0, 90]);
-
-	cube(10);
-	multmatrix(m_multiply(ma, mb))
-		translate([15, 0, 0])  cube(10);
-
-![m_multiply](images/lib-m_multiply-1.JPG)
-

docs/lib-m_translation.md

@@ -14,7 +14,7 @@ Generate a 4x4 transformation matrix which can pass into `multmatrix` to transla
 
 	cube(2, center = true); 
 	multmatrix(m_translation([5, 0, 0]))
-	    sphere(1,center = true);
+	    sphere(1);
 
 ![m_translation](images/lib-m_translation-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. You might have to set `twist` for connecting naturally. 
+- `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-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-torus_knot.md

@@ -10,7 +10,7 @@ Generate a path of [The (p,q)-torus knot](https://en.wikipedia.org/wiki/Torus_kn
 
 - `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. 
+- `phi_step` : The amount when increasing phi. 
 
 ## Examples
 
@@ -33,7 +33,7 @@ Generate a path of [The (p,q)-torus knot](https://en.wikipedia.org/wiki/Torus_kn
 		shape_pentagram_pts, 
 		concat(pts, [pts[0]]), 
 		closed = true,
-		twist = 188
+		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__/__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/along_with.scad

@@ -16,14 +16,12 @@ include <__private__/__to3d.scad>;
 // 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];
+
+    angles_defined = angles != undef;
 
     scale_step_vt = __is_float(scale) ? 
         scale_step() :
@@ -33,8 +31,20 @@ module along_with(points, angles, twist = 0, scale = 1.0) {
             scale[2] == undef ? 0 : (scale[2] - 1) / leng_points_minus_one
         ]; 
 
-    // get rotation matrice for sections
 
+    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],
@@ -42,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,
@@ -67,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]
         )
@@ -81,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) 
@@ -89,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)
@@ -98,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/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

@@ -8,12 +8,8 @@
 *
 **/
 
-function _bijection_edges_from(pts) = 
-    let(leng = len(pts))
-    concat(
-        [for(i = [0:leng - 2]) [pts[i], pts[i + 1]]], 
-        [[pts[len(pts) - 1], pts[0]]]
-    );
+include <__private__/__lines_from.scad>;
+include <__private__/__line_intersection.scad>;
     
 function _bijection_inward_edge_normal(edge) =  
     let(
@@ -45,30 +41,13 @@ function _bijection__bijection_offset_edges(edges, d) =
         )
         _bijection_offset_edge(edge, dx, dy)
     ];
-    
 
-function _bijection__bijection__bijection_offset_edges_intersection(edge1, edge2) = 
+function bijection_offset(pts, d, epsilon = 0.0001) = 
     let(
-        den = (edge2[1][1] - edge2[0][1]) * (edge1[1][0] - edge1[0][0]) - (edge2[1][0] - edge2[0][0]) * (edge1[1][1] - edge1[0][1])
-    )
-    // when den is 0, they are parallel or conincident edges
-    den == 0 ? [] : _bijection_offset__bijection__bijection__bijection_offset_edges_intersection_sub(edge1, edge2, den);
-
-function _bijection_offset__bijection__bijection__bijection_offset_edges_intersection_sub(edge1, edge2, den) = 
-    let(
-        ua = ((edge2[1][0] - edge2[0][0]) * (edge1[0][1] - edge2[0][1]) - (edge2[1][1] - edge2[0][1]) * (edge1[0][0] - edge2[0][0])) / den
-    ) 
-    [
-        edge1[0][0] + ua * (edge1[1][0] - edge1[0][0]),
-        edge1[0][1] + ua * (edge1[1][1] - edge1[0][1])
-    ];
-
-function bijection_offset(pts, d) = 
-    let(
-        es = _bijection_edges_from(pts), 
+        es = __lines_from(pts, true), 
         offset_es = _bijection__bijection_offset_edges(es, d),
         leng = len(offset_es),
-        last_p = _bijection__bijection__bijection_offset_edges_intersection(offset_es[leng - 1], offset_es[0])
+        last_p = __line_intersection(offset_es[leng - 1], offset_es[0], epsilon)
     )
     concat(
         [
@@ -76,7 +55,7 @@ function bijection_offset(pts, d) =
             let(
                 this_edge = offset_es[i],
                 next_edge = offset_es[i + 1],
-                p = _bijection__bijection__bijection_offset_edges_intersection(this_edge, next_edge)
+                p = __line_intersection(this_edge, next_edge, epsilon)
             )
             // p == p to avoid [nan, nan], because [nan, nan] != [nan, nan]
             if(p != [] && p == p) p

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(shapt_pts, true)
+    )
+    _in_shape_in_any_edges(edges, pt, epsilon) ? include_edge : 
+    _in_shape_sub(shapt_pts, leng, pt, false, leng - 1, 0);

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/parse_number.scad

@@ -4,14 +4,15 @@
 * @copyright Justin Lin, 2017
 * @license https://opensource.org/licenses/lgpl-3.0.html
 *
-* @see https://openhome.cc/eGossip/OpenSCAD/lib-parse_number.html
+* @see https://openhome.cc/eGossip/OpenSCAD/lib2-parse_number.html
 *
 **/ 
     
-function _str_to_int(t, i = 0, mapper = [["0", 0], ["1", 1], ["2", 2], ["3", 3], ["4", 4], ["5", 5], ["6", 6], ["7", 7], ["8", 8], ["9", 9]]) =  
-    i == len(mapper) ? -1 : (
-        mapper[i][0] == t ? mapper[i][1] : _str_to_int(t, i + 1)
-    );
+function _str_to_int(t) =  
+    let(
+        dict = [["0", 0], ["1", 1], ["2", 2], ["3", 3], ["4", 4], ["5", 5], ["6", 6], ["7", 7], ["8", 8], ["9", 9]],
+        n = dict[search(t, dict)[0]][1]
+    ) n;
     
 function _parse_positive_int(t, value = 0, i = 0) =
     i == len(t) ? value : _parse_positive_int(t, value * pow(10, i) + _str_to_int(t[i]), i + 1);

src/path_extrude.scad

@@ -16,151 +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_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
-        ];   
-
-    // 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/sub_str.scad

@@ -8,8 +8,9 @@
 *
 **/ 
 
-function sub_str(t, begin, end, result = "") =
-    end == undef ? sub_str(t, begin, len(t)) : (
-        begin == end ? result : sub_str(t, begin + 1, end, str(result, t[begin]))
-    );
+function _sub_str(t, begin, end) = 
+    begin == end ? "" : str(t[begin], sub_str(t, begin + 1, end));
+    
+function sub_str(t, begin, end) = 
+    end == undef ? _sub_str(t, begin, len(t)) : _sub_str(t, begin, end);
     

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