rewrite vnf_halfspace and fix bugs in half_of functions.

remove get_named_args

common.scad

@@ -407,7 +407,7 @@ function all_defined(v,recursive=false) =
 // Usage:
 //   anchr = get_anchor(anchor,center,<uncentered>,<dflt>);
 // Topics: Argument Handling
-// See Also: get_radius(), get_named_args()
+// See Also: get_radius()
 // Description:
 //   Calculated the correct anchor from `anchor` and `center`.  In order:
 //   - If `center` is not `undef` and `center` evaluates as true, then `CENTER` (`[0,0,0]`) is returned.
@@ -437,7 +437,7 @@ function get_anchor(anchor,center,uncentered=BOT,dflt=CENTER) =
 // Usage:
 //   r = get_radius(<r1=>, <r2=>, <r=>, <d1=>, <d2=>, <d=>, <dflt=>);
 // Topics: Argument Handling
-// See Also: get_anchor(), get_named_args()
+// See Also: get_anchor()
 // Description:
 //   Given various radii and diameters, returns the most specific radius.  If a diameter is most
 //   specific, returns half its value, giving the radius.  If no radii or diameters are defined,
@@ -480,88 +480,6 @@ function get_radius(r1, r2, r, d1, d2, d, dflt) =
     : dflt;
 
 
-// Function: get_named_args()
-// Usage:
-//   function f(pos1=_UNDEF, pos2=_UNDEF,...,named1=_UNDEF, named2=_UNDEF, ...) = let(args = get_named_args([pos1, pos2, ...], [[named1, default1], [named2, default2], ...]), named1=args[0], named2=args[1], ...)
-// Topics: Argument Handling
-// See Also: get_anchor(), get_radius()
-// Description:
-//   Given the values of some positional and named arguments, returns a list of the values assigned to
-//   named parameters.  in the following steps:
-//   - First, all named parameters which were explicitly assigned in the function call take their
-//     provided value.
-//   - Then, any positional arguments are assigned to remaining unassigned
-//     parameters; this is governed both by the `priority` entries (if there are `N` positional
-//     arguments, then the `N` parameters with lowest `priority` value will be assigned) and by the
-//     order of the positional arguments (matching that of the assigned named parameters).  If no
-//     priority is given, then these two ordering coincide: parameters are assigned in order, starting
-//     from the first one.
-//   - Finally, any remaining named parameters can take default values.  If no default values are
-//     given, then `undef` is used.
-//   .
-//   This allows an author to declare a function prototype with named or optional parameters, so that
-//   the user may then call this function using either positional or named parameters. In practice the
-//   author will declare the function as using *both* positional and named parameters, and let
-//   `get_named_args()` do the parsing from the whole set of arguments.  See the example below.
-//   .
-//   This supports the user explicitly passing `undef` as a function argument.  To distinguish between
-//   an intentional `undef` and the absence of an argument, we use a custom `_UNDEF` value as a guard
-//   marking the absence of any arguments (in practice, `_UNDEF` is a random-generated string, which
-//   will never coincide with any useful user value).  This forces the author to declare all the
-//   function parameters as having `_UNDEF` as their default value.
-// Arguments:
-//   positional = The list of values of positional arguments.
-//   named = The list of named arguments; each entry of the list has the form `[passed-value, <default-value>, <priority>]`, where `passed-value` is the value that was passed at function call; `default-value` is the value that will be used if nothing is read from either named or positional arguments; `priority` is the priority assigned to this argument (lower means more priority, default value is `+inf`). Since stable sorting is used, if no priority at all is given, all arguments will be read in order.
-//   _undef = The default value used by the calling function for all arguments. The default value, `_UNDEF`, is a random string. This value **must** be the default value of all parameters in the outer function call (see example below).
-//
-// Example: a function with prototype `f(named1,< <named2>, named3 >)`
-//   function f(_p1=_UNDEF, _p2=_UNDEF, _p3=_UNDEF,
-//              arg1=_UNDEF, arg2=_UNDEF, arg3=_UNDEF) =
-//      let(named = get_named_args([_p1, _p2, _p3],
-//          [[arg1, "default1",0], [arg2, "default2",2], [arg3, "default3",1]]))
-//      named;
-//   // all default values or all parameters provided:
-//   echo(f());
-//   // ["default1", "default2", "default3"]
-//   echo(f("given2", "given3", arg1="given1"));
-//   // ["given1", "given2", "given3"]
-//   
-//   // arg1 has highest priority, and arg3 is higher than arg2:
-//   echo(f("given1"));
-//   // ["given1", "default2", "default3"]
-//   echo(f("given3", arg1="given1"));
-//   // ["given1", "default2", "given3"]
-//   
-//   // explicitly passing undef is allowed:
-//   echo(f(undef, arg1="given1", undef));
-//   // ["given1", undef, undef]
-
-/* Note: however tempting it might be, it is *not* possible to accept
- * named argument as a list [named1, named2, ...] (without default
- * values), because the values [named1, named2...] themselves might be
- * lists, and we will not be able to distinguish the two cases. */
-function get_named_args(positional, named, _undef=_UNDEF) =
-    let(deft = [for(p=named) p[1]], // default is undef
-        // indices of the values to fetch from positional args:
-        unknown = [for(x=enumerate(named)) if(x[1][0]==_undef) x[0]],
-        // number of values given to positional arguments:
-        n_positional = count_true([for(p=positional) p!=_undef]))
-    assert(n_positional <= len(unknown),
-      str("too many positional arguments (", n_positional, " given, ",
-          len(unknown), " required)"))
-    let(
-        // those elements which have no priority assigned go last (prio=+∞):
-        prio = sortidx([for(u=unknown) default(named[u][2], 1/0)]),
-        // list of indices of values assigned from positional arguments:
-        assigned = [for(a=sort([for(i=[0:1:n_positional-1]) prio[i]]))
-          unknown[a]])
-    [ for(e = enumerate(named))
-      let(idx=e[0], val=e[1][0], ass=search(idx, assigned))
-        val != _undef ? val :
-        ass != [] ? positional[ass[0]] :
-        deft[idx] ];
-
-
 // Function: scalar_vec3()
 // Usage:
 //   vec = scalar_vec3(v, <dflt>);

mutators.scad

@@ -96,23 +96,23 @@ module bounding_box(excess=0, planar=false) {
 // Function&Module: half_of()
 //
 // Usage: as module
-//   half_of(v, <cp>, <s>) ...
+//   half_of(v, <cp>, <s>, <planar>) ...
 // Usage: as function
-//   half_of(v, <cp>, p, <s>)...
+//   result = half_of(p,v,<cp>);
 //
 // Description:
-//   Slices an object at a cut plane, and masks away everything that is on one side.
+//   Slices an object at a cut plane, and masks away everything that is on one side.  The v parameter is either a plane specification or
-//   * Called as a function with a path in the `p` argument, returns the
+//   a normal vector.  The s parameter is needed for the module
-//       intersection of path `p` and given half-space.
+//   version to control the size of the masking cube, which affects preview display.
-//   * Called as a function with a 2D path in the `p` argument
+//   When called as a function, you must supply a vnf, path or region in p.  If planar is set to true for the module version the operation
-//       and a 2D vector `p`, returns the intersection of path `p` and given
+//   is performed in  and UP and DOWN are treated as equivalent to BACK and FWD respectively.  
-//       half-plane.
 //
 // Arguments:
+//   p = path, region or VNF to slice.  (Function version)
 //   v = Normal of plane to slice at.  Keeps everything on the side the normal points to.  Default: [0,0,1] (UP)
-//   cp = If given as a scalar, moves the cut plane along the normal by the given amount.  If given as a point, specifies a point on the cut plane.  This can be used to shift where it slices the object at.  Default: [0,0,0]
+//   cp = If given as a scalar, moves the cut plane along the normal by the given amount.  If given as a point, specifies a point on the cut plane.  Default: [0,0,0]
-//   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, it messes with centering your view.  Default: 100
+//   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, it messes with centering your view.   Ignored for function version.   Default: 1000
-//   planar = If true, this becomes a 2D operation.  When planar, a `v` of `UP` or `DOWN` becomes equivalent of `BACK` and `FWD` respectively.
+//   planar = If true, perform a 2D operation.  When planar, a `v` of `UP` or `DOWN` becomes equivalent of `BACK` and `FWD` respectively.
 //
 // Examples:
 //   half_of(DOWN+BACK, cp=[0,-10,0]) cylinder(h=40, r1=10, r2=0, center=false);
@@ -147,23 +147,43 @@ module half_of(v=UP, cp, s=1000, planar=false)
     }
 }
 
-function half_of(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF, _arg4=_UNDEF,
+function half_of(p, v=UP, cp) =
-    v=_UNDEF, cp=_UNDEF, p=_UNDEF, s=_UNDEF) =
+    is_vnf(p) ?
-    let(args=get_named_args([_arg1, _arg2, _arg3, _arg4],
+       assert(is_vector(v) && (len(v)==3 || len(v)==4),str("Must give 3-vector or plane specification",v))
-        [[v,undef,0], [cp,0,2], [p,undef,1], [s, 1e4]]),
+       assert(select(v,0,2)!=[0,0,0], "vector v must be nonzero")
-        v=args[0], cp0=args[1], p=args[2], s=args[3],
+       let(
-        cp = is_num(cp0) ? cp0*unit(v) : cp0)
+            plane = is_vector(v,4) ? assert(cp==undef, "Don't use cp with plane definition.") v
-    assert(is_vector(v,2)||is_vector(v,3),
+                  : is_undef(cp) ? [each v, 0]
-      "must provide a half-plane or half-space")
+                  : is_num(cp) ? [each v, cp*(v*v)/norm(v)]
-    let(d=len(v))
+                  : assert(is_vector(cp,3),"Centerpoint must be a 3-vector")
-    assert(len(cp) == d, str("cp must have dimension ", d))
+                    [each v, cp*v]
-    is_vector(p) ?
+       )
-        assert(len(p) == d, str("vector must have dimension ", d))
+       vnf_halfspace(plane, p)
-        let(z=(p-cp)*v) (z >= 0 ? p : p - (z*v)/(v*v))
+   : is_path(p) || is_region(p) ?
-        :
+      let(
-    p == [] ? [] : // special case: empty path remains empty
+          v = (v==UP)? BACK : (v==DOWN)? FWD : v,
+          cp = is_undef(cp) ? [0,0]
+             : is_num(cp) ? v*cp 
+             : assert(is_vector(cp,2) || (is_vector(cp,3) && cp.z==0),"Centerpoint must be 2-vector")
+               cp
+      )
+      assert(is_vector(v,2) || (is_vector(v,3) && v.z==0),"Must give 2-vector")
+      assert(!all_zero(v), "Vector v must be nonzero")
+      let(
+          bounds = pointlist_bounds(move(-cp,p)),
+          L = 2*max(flatten(bounds)),
+          n = unit(v),
+          u = [-n.y,n.x],
+          box = [cp+u*L, cp+(v+u)*L, cp+(v-u)*L, cp-u*L]
+      )
+      intersection(box,p)
+   : assert(false, "Input must be a region, path or VNF");
+
+          
+
+/*  This code cut 3d paths but leaves behind connecting line segments
     is_path(p) ?
-        assert(len(p[0]) == d, str("path must have dimension ", d))
+        //assert(len(p[0]) == d, str("path must have dimension ", d))
         let(z = [for(x=p) (x-cp)*v])
         [ for(i=[0:len(p)-1]) each concat(z[i] >= 0 ? [p[i]] : [],
             // we assume a closed path here;
@@ -174,17 +194,8 @@ function half_of(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF, _arg4=_UNDEF,
             // create self-intersection or whiskers:
             z[i]*z[j] >= 0 ? [] : [(z[j]*p[i]-z[i]*p[j])/(z[j]-z[i])]) ]
         :
-    is_vnf(p) ?
+*/
-    // we must put is_vnf() before is_region(), because most triangulated
+
-    // VNFs will pass is_region() test
-    vnf_halfspace(halfspace=concat(v,[-v*cp]), vnf=p) :
-    is_region(p) ?
-        assert(len(v) == 2, str("3D vector not compatible with region"))
-        let(u=unit(v), w=[-u[1], u[0]],
-            R=[[cp+s*w, cp+s*(v+v), cp+s*(v-w), cp-s*w]]) // half-plane
-        intersection(R, p)
-        :
-    assert(false, "must pass either a point, a path, a region, or a VNF");
 
 // Function&Module: left_half()
 //
@@ -192,17 +203,16 @@ function half_of(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF, _arg4=_UNDEF,
 //   left_half(<s>, <x>) ...
 //   left_half(planar=true, <s>, <x>) ...
 // Usage: as function
-//   left_half(<s>, <x>, path)
+//   result = left_half(p, <x>);
-//   left_half(<s>, <x>, region)
-//   left_half(<s>, <x>, vnf)
 //
 // Description:
 //   Slices an object at a vertical Y-Z cut plane, and masks away everything that is right of it.
 //
 // Arguments:
+//   p = VNF, region or path to slice (function version)
 //   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, OpenSCAD's preview rendering may be incorrect.  Default: 10000
 //   x = The X coordinate of the cut-plane.  Default: 0
-//   planar = If true, this becomes a 2D operation.
+//   planar = If true, perform a 2D operation.
 //
 // Examples:
 //   left_half() sphere(r=20);
@@ -223,28 +233,26 @@ module left_half(s=1000, x=0, planar=false)
         }
     }
 }
-function left_half(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF,
+function left_half(p,x=0) = half_of(p, LEFT, [x,0,0]);
-    x=_UNDEF, p=_UNDEF, s=_UNDEF) =
-    let(args=get_named_args([_arg1, _arg2, _arg3],
-    [[x, 0,1], [p,undef,0], [s, 1e4]]),
-        x=args[0], p=args[1], s=args[2])
-    half_of(v=[-1,0,0], cp=x, p=p);
 
 
 
 // Function&Module: right_half()
 //
-// Usage:
+// Usage: as module
 //   right_half([s], [x]) ...
 //   right_half(planar=true, [s], [x]) ...
+// Usage: as function
+//   result = right_half(p, <x>);
 //
 // Description:
 //   Slices an object at a vertical Y-Z cut plane, and masks away everything that is left of it.
 //
 // Arguments:
+//   p = VNF, region or path to slice (function version)
 //   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, OpenSCAD's preview rendering may be incorrect.  Default: 10000
 //   x = The X coordinate of the cut-plane.  Default: 0
-//   planar = If true, this becomes a 2D operation.
+//   planar = If true perform a 2D operation.
 //
 // Examples(FlatSpin,VPD=175):
 //   right_half() sphere(r=20);
@@ -265,28 +273,26 @@ module right_half(s=1000, x=0, planar=false)
         }
     }
 }
-function right_half(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF,
+function right_half(p,x=0) = half_of(p, RIGHT, [x,0,0]);
-    x=_UNDEF, p=_UNDEF, s=_UNDEF) =
-    let(args=get_named_args([_arg1, _arg2, _arg3],
-    [[x, 0,1], [p,undef,0], [s, 1e4]]),
-        x=args[0], p=args[1], s=args[2])
-    half_of(v=[1,0,0], cp=x, p=p);
 
 
 
 // Function&Module: front_half()
 //
 // Usage:
-//   front_half([s], [y]) ...
+//   front_half(<s>, <y>) ...
-//   front_half(planar=true, [s], [y]) ...
+//   front_half(planar=true, <s>, <y>) ...
+// Usage: as function
+//   result = front_half(p, <y>);
 //
 // Description:
 //   Slices an object at a vertical X-Z cut plane, and masks away everything that is behind it.
 //
 // Arguments:
+//   p = VNF, region or path to slice (function version)
 //   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, OpenSCAD's preview rendering may be incorrect.  Default: 10000
 //   y = The Y coordinate of the cut-plane.  Default: 0
-//   planar = If true, this becomes a 2D operation.
+//   planar = If true perform a 2D operation.
 //
 // Examples(FlatSpin,VPD=175):
 //   front_half() sphere(r=20);
@@ -307,28 +313,26 @@ module front_half(s=1000, y=0, planar=false)
         }
     }
 }
-function front_half(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF,
+function front_half(p,y=0) = half_of(p, FRONT, [0,y,0]);
-    x=_UNDEF, p=_UNDEF, s=_UNDEF) =
-    let(args=get_named_args([_arg1, _arg2, _arg3],
-    [[x, 0,1], [p,undef,0], [s, 1e4]]),
-        x=args[0], p=args[1], s=args[2])
-    half_of(v=[0,-1,0], cp=x, p=p);
 
 
 
 // Function&Module: back_half()
 //
 // Usage:
-//   back_half([s], [y]) ...
+//   back_half(<s>, <y>) ...
-//   back_half(planar=true, [s], [y]) ...
+//   back_half(planar=true, <s>, <y>) ...
+// Usage: as function
+//   result = back_half(p, <y>);
 //
 // Description:
 //   Slices an object at a vertical X-Z cut plane, and masks away everything that is in front of it.
 //
 // Arguments:
+//   p = VNF, region or path to slice (function version)
 //   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, OpenSCAD's preview rendering may be incorrect.  Default: 10000
 //   y = The Y coordinate of the cut-plane.  Default: 0
-//   planar = If true, this becomes a 2D operation.
+//   planar = If true perform a 2D operation.
 //
 // Examples:
 //   back_half() sphere(r=20);
@@ -349,24 +353,22 @@ module back_half(s=1000, y=0, planar=false)
         }
     }
 }
-function back_half(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF,
+function back_half(p,y=0) = half_of(p, BACK, [0,y,0]);
-    x=_UNDEF, p=_UNDEF, s=_UNDEF) =
-    let(args=get_named_args([_arg1, _arg2, _arg3],
-    [[x, 0,1], [p,undef,0], [s, 1e4]]),
-        x=args[0], p=args[1], s=args[2])
-    half_of(v=[0,1,0], cp=x, p=p);
 
 
 
 // Function&Module: bottom_half()
 //
 // Usage:
-//   bottom_half([s], [z]) ...
+//   bottom_half(<s>, <z>) ...
+// Usage: as function
+//   result = bottom_half(p, <z>);
 //
 // Description:
 //   Slices an object at a horizontal X-Y cut plane, and masks away everything that is above it.
 //
 // Arguments:
+//   p = VNF, region or path to slice (function version)
 //   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, OpenSCAD's preview rendering may be incorrect.  Default: 10000
 //   z = The Z coordinate of the cut-plane.  Default: 0
 //
@@ -383,24 +385,21 @@ module bottom_half(s=1000, z=0)
         }
     }
 }
-function bottom_half(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF,
+function bottom_half(p,z=0) = half_of(p,BOTTOM,[0,0,z]);
-    x=_UNDEF, p=_UNDEF, s=_UNDEF) =
-    let(args=get_named_args([_arg1, _arg2, _arg3],
-    [[x, 0,1], [p,undef,0], [s, 1e4]]),
-        x=args[0], p=args[1], s=args[2])
-    half_of(v=[0,0,-1], cp=x, p=p);
 
 
 
 // Function&Module: top_half()
 //
 // Usage:
-//   top_half([s], [z]) ...
+//   top_half(<s>, <z>) ...
+//   result = top_half(p, <z>);
 //
 // Description:
 //   Slices an object at a horizontal X-Y cut plane, and masks away everything that is below it.
 //
 // Arguments:
+//   p = VNF, region or path to slice (function version)
 //   s = Mask size to use.  Use a number larger than twice your object's largest axis.  If you make this too large, OpenSCAD's preview rendering may be incorrect.  Default: 10000
 //   z = The Z coordinate of the cut-plane.  Default: 0
 //
@@ -417,12 +416,7 @@ module top_half(s=1000, z=0)
         }
     }
 }
-function top_half(_arg1=_UNDEF, _arg2=_UNDEF, _arg3=_UNDEF,
+function top_half(p,z=0) = half_of(p,UP,[0,0,z]);
-    x=_UNDEF, p=_UNDEF, s=_UNDEF) =
-    let(args=get_named_args([_arg1, _arg2, _arg3],
-    [[x, 0,1], [p,undef,0], [s, 1e4]]),
-        x=args[0], p=args[1], s=args[2])
-    half_of(v=[0,0,1], cp=x, p=p);
 
 
 

vnf.scad

@@ -1031,144 +1031,111 @@ module vnf_validate(vnf, size=1, show_warns=true, check_isects=false) {
 
 // Function: vnf_halfspace()
 // Usage:
-//   vnf_halfspace([a,b,c,d], vnf)
+//   newvnf = vnf_halfspace(plane, vnf, <closed>);
 // Description:
-//   returns the intersection of the VNF with the given half-space.
+//   Returns the intersection of the vnf with a half space.  The half space is defined by
+//   plane = [A,B,C,D], taking the side where the normal [A,B,C] points: Ax+By+Cz≥D.
+//   If closed is set to false then the cut face is not included in the vnf.  This could
+//   allow further extension of the vnf by merging with other vnfs.  
 // Arguments:
-//   halfspace = half-space to intersect with, given as the four coefficients of the affine inequation a\*x+b\*y+c\*z≥ d.
+//   plane = plane defining the boundary of the half space
+//   vnf = vnf to cut
+//   closed = if false do not return include cut face(s).  Default: true
+// Example:
+//   vnf = cube(10,center=true);
+//   cutvnf = vnf_halfspace([-1,1,-1,0], vnf);
+//   vnf_polyhedron(cutvnf);
+// Example:  Cut face has 2 components
+//   vnf = path_sweep(circle(r=4, $fn=16),
+//                    circle(r=20, $fn=64),closed=true);
+//   cutvnf = vnf_halfspace([-1,1,-4,0], vnf);
+//   vnf_polyhedron(cutvnf);
+// Example: Cut face is not simply connected
+//   vnf = path_sweep(circle(r=4, $fn=16),
+//                    circle(r=20, $fn=64),closed=true);
+//   cutvnf = vnf_halfspace([0,0.7,-4,0], vnf);
+//   vnf_polyhedron(cutvnf);
+// Example: Cut object has multiple components
+//   function knot(a,b,t) =   // rolling knot 
+//        [ a * cos (3 * t) / (1 - b* sin (2 *t)), 
+//          a * sin( 3 * t) / (1 - b* sin (2 *t)), 
+//        1.8 * b * cos (2 * t) /(1 - b* sin (2 *t))]; 
+//   a = 0.8; b = sqrt (1 - a * a); 
+//   ksteps = 400;
+//   knot_path = [for (i=[0:ksteps-1]) 50 * knot(a,b,(i/ksteps)*360)];
+//   ushape = [[-10, 0],[-10, 10],[ -7, 10],[ -7, 2],[  7, 2],[  7, 7],[ 10, 7],[ 10, 0]];
+//   knot=path_sweep(ushape, knot_path, closed=true, method="incremental");
+//   cut_knot = vnf_halfspace([1,0,0,0], knot);
+//   vnf_polyhedron(cut_knot);
+function vnf_halfspace(plane, vnf, closed=true) =
+    let(
+         inside = [for(x=vnf[0]) plane*[each x,-1] >= 0 ? 1 : 0],
+         vertexmap = [0,each cumsum(inside)],
+         faces_edges_vertices = _vnfcut(plane, vnf[0],vertexmap,inside, vnf[1], last(vertexmap)),
+         newvert = concat(bselect(vnf[0],inside), faces_edges_vertices[2])
+    )
+    closed==false ? [newvert, faces_edges_vertices[0]] :
+    let(
+        allpaths = _assemble_paths(newvert, faces_edges_vertices[1]),
+        newpaths = [for(p=allpaths) if (len(p)>=3) p
+                                    else assert(approx(p[0],p[1]),"Orphan edge found when assembling cut edges.")
+           ]
+    )
+    len(newpaths)<=1 ? [newvert, concat(faces_edges_vertices[0], newpaths)] 
+    :
+      let(
+           faceregion = [for(p=newpaths) project_plane(select(newvert,p), plane)],
+           facevnf = region_faces(faceregion,reverse=true),
+           faceverts = lift_plane(facevnf[0], plane)
+      )
+      vnf_merge([[newvert, faces_edges_vertices[0]], [faceverts, facevnf[1]]]);
 
-function _vnf_halfspace_pts(halfspace, points, faces,
+
-  inside=undef, coords=[], map=[]) =
+function _assemble_paths(vertices, edges, paths=[],i=0) =
-/* Recursive function to compute the intersection of points (and edges,
+     i==len(edges) ? paths :
- * but not faces) with with the half-space.
+     norm(vertices[edges[i][0]]-vertices[edges[i][1]])<EPSILON ? echo(degen=i)_assemble_paths(vertices,edges,paths,i+1) :
- * Parameters:
+     let(    // Find paths that connects on left side and right side of the edges (if one exists)
- * halfspace  a vector(4)
+         left = [for(j=idx(paths)) if (approx(vertices[last(paths[j])],vertices[edges[i][0]])) j],
- * points     a list of points3d
+         right = [for(j=idx(paths)) if (approx(vertices[edges[i][1]],vertices[paths[j][0]])) j]
- * faces      a list of indexes in points
+     )
- * inside     a vector{bool} determining which points belong to the
+     assert(len(left)<=1 && len(right)<=1)
- *            half-space; if undef, it is initialized at first loop.
+     let(              
- * coords     the coordinates of the points in the intersection
+          keep_path = list_remove(paths,concat(left,right)),
- * map        the logical map (old point) → (new point(s)):
+          update_path = left==[] && right==[] ? edges[i] 
- *   if point i is kept, then map[i] = new-index-for-i;
+                      : left==[] ? concat([edges[i][0]],paths[right[0]])
- *   if point i is dropped, then map[i] = [[j1, k1], [j2, k2], …],
+                      : right==[] ? concat(paths[left[0]],[edges[i][1]])
- *      where points j1,… are kept (old index)
+                      : left != right ? concat(paths[left[0]], paths[right[0]])
- *      and k1,… are the matching intersections (new index).
+                      : paths[left[0]]
- * Returns the triple [coords, map, inside].
+     )
- *
+     _assemble_paths(vertices, edges, concat(keep_path, [update_path]), i+1);
- */
+
-    let(i=len(map), n=len(coords)) // we are currently processing point i
+
-    // termination test:
+function _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount, newfaces=[], newedges=[], newvertices=[], i=0) =
-    i >= len(points) ? [ coords, map, inside ] :
+   i==len(faces) ? [newfaces, newedges, newvertices] :
-    let(inside = !is_undef(inside) ? inside :
+   let(
-        [for(x=points) halfspace*concat(x,[-1]) >= 0],
+        pts_inside = select(inside,faces[i])
-        pi = points[i])
+   )
-    // inside half-space: keep the point (and reindex)
+   all(pts_inside) ? _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount,
-    inside[i] ? _vnf_halfspace_pts(halfspace, points, faces, inside,
+                             concat(newfaces, [select(vertexmap,faces[i])]), newedges, newvertices, i+1):
-        concat(coords, [pi]), concat(map, [n]))
+   !any(pts_inside) ? _vnfcut(plane, vertices, vertexmap,inside, faces, vertcount, newfaces, newedges, newvertices, i+1):
-    : // else: compute adjacent vertices (adj)
+   let(
-    let(adj = unique([for(f=faces) let(m=len(f), j=search(i, f)[0])
+        first = search([[1,0]],pair(pts_inside,wrap=true),0)[0],
-      each if(j!=undef) [f[(j+1)%m], f[(j+m-1)%m]] ]),
+        second = search([[0,1]],pair(pts_inside,wrap=true),0)[0]
-    // filter those which lie in half-space:
+   )
-        adj2 = [for(x=adj) if(inside[x]) x],
+   assert(len(first)==1 && len(second)==1, "Found concave face in VNF.  Run vnf_triangulate first to ensure convex faces.")
-        zi = halfspace*concat(pi, [-1]))
+   let(
-    _vnf_halfspace_pts(halfspace, points, faces, inside,
+        newface = [each select(vertexmap,select(faces[i],second[0]+1,first[0])),vertcount, vertcount+1],
-        // new points: we append all these intersection points
+        newvert = [plane_line_intersection(plane, select(vertices,select(faces[i],first[0],first[0]+1)),eps=0),
-        concat(coords, [for(j=adj2) let(zj=halfspace*concat(points[j],[-1]))
+                   plane_line_intersection(plane, select(vertices,select(faces[i],second[0],second[0]+1)),eps=0)]
-            (zi*points[j]-zj*pi)/(zi-zj)]),
+   )
-        // map: we add the info
+   true //!approx(newvert[0],newvert[1])
-        concat(map, [[for(y=enumerate(adj2)) [y[1], n+y[0]]]]));
+       ? _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount+2,
-function _vnf_halfspace_face(face, map, inside, i=0,
+                 concat(newfaces, [newface]), concat(newedges,[[vertcount+1,vertcount]]),concat(newvertices,newvert),i+1)
-    newface=[], newedge=[], exit) =
+   :len(newface)>3
-/* Recursive function to intersect a face of the VNF with the half-plane.
+       ? _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount+1,
- * Arguments:
+                 concat(newfaces, [list_head(newface)]), newedges,concat(newvertices,[newvert[0]]),i+1)
- *   face: the list of points of the face (old indices).
+   :
- *   map: as produced by _vnf_halfspace_pts
+   _vnfcut(plane, vertices, vertexmap, inside, faces, vertcount,newfaces, newedges, newvert, i+1);
- *   inside: vector{bool} containing half-space info
+ 
- *   i: index for iteration
- *   exit: boolean; is first point in newedge an exit or an entrance from
- *     half-space?
- *   newface: list of (new indexes of) points on the face
- *   newedge: list of new points on the plane (even number of points)
- *  Return value: [newface, new-edges], where new-edges is a list of
- *  pairs [entrance-node, exit-node] (new indices).
- */
-// termination condition:
-    (i >= len(face)) ? [ newface,
-    // if exit==true then we return newedge[1,0], newedge[3,2], ...
-    // otherwise newedge[0,1], newedge[2,3], ...;
-    // all edges are oriented (entrance->exit), so that by following the
-    // arrows we obtain a correctly-oriented face:
-    let(k = exit ? 0 : 1)
-    [for(i=[0:2:len(newedge)-2]) [newedge[i+k], newedge[i+1-k]]] ]
-    : // recursion case: p is current point on face, q is next point
-    let(p = face[i], q = face[(i+1)%len(face)],
-        // if p is inside half-plane, keep it in the new face:
-        newface0 = inside[p] ?  concat(newface, [map[p]]) : newface)
-        // if the current segment does not intersect, this is all:
-        inside[p] == inside[q] ? _vnf_halfspace_face(face, map, inside, i+1,
-            newface0, newedge, exit)
-        : // otherwise, we must add the intersection point:
-        // rename the two points p,q as inner and outer point:
-        let(in = inside[p] ? p : q, out = p+q-in,
-            inter=[for(a=map[out]) if(a[0]==in) a[1]][0])
-        _vnf_halfspace_face(face, map, inside, i+1,
-            concat(newface0, [inter]),
-            concat(newedge, [inter]),
-            is_undef(exit) ? inside[p] : exit);
-function _vnf_halfspace_path_search_edge(edge, paths, i=0, ret=[undef,undef]) =
-/* given an oriented edge [x,y] and a set of oriented paths,
- * returns the indices [i,j] of paths [before, after] given edge
- */
-    // termination condition
-    i >= len(paths) ? ret:
-    _vnf_halfspace_path_search_edge(edge, paths, i+1,
-       [last(paths[i]) == edge[0] ? i : ret[0],
-        paths[i][0] == edge[1] ? i : ret[1]]);
-function _vnf_halfspace_paths(edges, i=0, paths=[]) =
-/* given a set of oriented edges [x,y],
-   returns all paths [x,y,z,..] that may be formed from these edges.
-   A closed path will be returned with equal first and last point.
-   i: index of currently examined edge
- */
-    i >= len(edges) ? paths : // termination condition
-    let(e=edges[i], s = _vnf_halfspace_path_search_edge(e, paths))
-        _vnf_halfspace_paths(edges, i+1,
-        // we keep all paths untouched by e[i]
-        concat([for(i=[0:1:len(paths)-1]) if(i!= s[0] && i != s[1]) paths[i]],
-        is_undef(s[0])? (
-            // fresh e: create a new path
-            is_undef(s[1]) ? [e] :
-            // e attaches to beginning of previous path
-            [concat([e[0]], paths[s[1]])]
-        ) :// edge attaches to end of previous path
-        is_undef(s[1]) ? [concat(paths[s[0]], [e[1]])] :
-        // edge merges two paths
-        s[0] != s[1] ? [concat(paths[s[0]], paths[s[1]])] :
-        // edge closes a loop
-        [concat(paths[s[0]], [e[1]])]));
-function vnf_halfspace(_arg1=_UNDEF, _arg2=_UNDEF,
-    halfspace=_UNDEF, vnf=_UNDEF) =
-    // here is where we wish that OpenSCAD had array lvalues...
-    let(args=get_named_args([_arg1, _arg2], [[halfspace],[vnf]]),
-        halfspace=args[0], vnf=args[1])
-    assert(is_vector(halfspace, 4),
-        "half-space must be passed as a length 4 affine form")
-    assert(is_vnf(vnf), "must pass a vnf")
-        // read points
-    let(tmp1=_vnf_halfspace_pts(halfspace, vnf[0], vnf[1]),
-        coords=tmp1[0], map=tmp1[1], inside=tmp1[2],
-        // cut faces and generate edges
-        tmp2= [for(f=vnf[1]) _vnf_halfspace_face(f, map, inside)],
-        newfaces=[for(x=tmp2) if(x[0]!=[]) x[0]],
-        newedges=[for(x=tmp2) each x[1]],
-        // generate new faces
-        paths=_vnf_halfspace_paths(newedges),
-        reg = [for(p=paths) project_plane(select(coords,p), halfspace)],
-        regvnf = region_faces(reg,reverse=true),
-        regvert = lift_plane(regvnf[0], halfspace)
-        //loops=[for(p=paths) if(coords[p[0]] == coords[last(p)]) reverse(p)])
-        )
-        vnf_merge([[coords, newfaces], [regvert, regvnf[1]]]);
-//    [coords, concat(newfaces, loops)];
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap