1//////////////////////////////////////////////////////////////////////
   2// LibFile: beziers.scad
   3//   Bezier curves and surfaces are way to represent smooth curves and smoothly curving
   4//   surfaces with a set of control points.  The curve or surface is defined by
   5//   the control points, but usually only passes through the first and last control point (the endpoints).
   6//   This file provides some
   7//   aids to constructing the control points, and highly optimized functions for
   8//   computing the Bezier curves and surfaces given by the control points, 
   9// Includes:
  10//   include <BOSL2/std.scad>
  11//   include <BOSL2/beziers.scad>
  12// FileGroup: Advanced Modeling
  13// FileSummary: Bezier curves and surfaces.
  14//////////////////////////////////////////////////////////////////////
  15
  16// Terminology:
  17//   Path = A series of points joined by straight line segements.
  18//   Bezier Curve = A polynomial curve defined by a list of control points.  The curve starts at the first control point and ends at the last one.  The other control points define the shape of the curve and they are often *NOT* on the curve
  19//   Control Point = A point that influences the shape of the Bezier curve.
  20//   Degree = The degree of the polynomial used to make the bezier curve.  A bezier curve of degree N will have N+1 control points.  Most beziers are cubic (degree 3).  The higher the degree, the more the curve can wiggle.  
  21//   Bezier Parameter = A parameter, usually `u` below, that ranges from 0 to 1 to trace out the bezier curve.  When `u=0` you get the first control point and when `u=1` you get the last control point. Intermediate points are traced out *non-uniformly*.  
  22//   Bezier Path = A list of bezier control points corresponding to a series of Bezier curves that connect together, end to end.  Because they connect, the endpoints are shared between control points and are not repeated, so a degree 3 bezier path representing two bezier curves will have seven entries to represent two sets of four control points.    **NOTE:** A "bezier path" is *NOT* a standard path
  23//   Bezier Patch = A two-dimensional arrangement of Bezier control points that generate a bounded curved Bezier surface.  A Bezier patch is a (N+1) by (M+1) grid of control points, which defines surface with four edges (in the non-degenerate case). 
  24//   Bezier Surface = A surface defined by a list of one or more bezier patches.
  25//   Spline Steps = The number of straight-line segments used to approximate a Bezier curve.  The more spline steps, the better the approximation to the curve, but the slower it will be to generate.  This plays a role analogous to `$fn` for circles.  Usually defaults to 16.
  26
  27
  28// Section: Bezier Curves
  29
  30// Function: bezier_points()
  31// Synopsis: Computes one or more specified points along a bezier curve.
  32// SynTags: Path
  33// Topics: Bezier Curves
  34// See Also: bezier_curve(), bezier_curvature(), bezier_tangent(), bezier_derivative(), bezier_points()
  35// Usage:
  36//   pt = bezier_points(bezier, u);
  37//   ptlist = bezier_points(bezier, RANGE);
  38//   ptlist = bezier_points(bezier, LIST);
  39// Description:
  40//   Computes points on a bezier curve with control points specified by `bezier` at parameter values
  41//   specified by `u`, which can be a scalar or a list.  The value `u=0` gives the first endpoint; `u=1` gives the final endpoint,
  42//   and intermediate values of `u` fill in the curve in a non-uniform fashion.  This function uses an optimized method which
  43//   is best when `u` is a long list and the bezier degree is 10 or less.  The degree of the bezier
  44//   curve is `len(bezier)-1`.
  45// Arguments:
  46//   bezier = The list of endpoints and control points for this bezier curve.
  47//   u = Parameter values for evaluating the curve, given as a single value, a list or a range.  
  48// Example(2D): Quadratic (Degree 2) Bezier.
  49//   bez = [[0,0], [30,30], [80,0]];
  50//   debug_bezier(bez, N=len(bez)-1);
  51//   translate(bezier_points(bez, 0.3)) color("red") sphere(1);
  52// Example(2D): Cubic (Degree 3) Bezier
  53//   bez = [[0,0], [5,35], [60,-25], [80,0]];
  54//   debug_bezier(bez, N=len(bez)-1);
  55//   translate(bezier_points(bez, 0.4)) color("red") sphere(1);
  56// Example(2D): Degree 4 Bezier.
  57//   bez = [[0,0], [5,15], [40,20], [60,-15], [80,0]];
  58//   debug_bezier(bez, N=len(bez)-1);
  59//   translate(bezier_points(bez, 0.8)) color("red") sphere(1);
  60// Example(2D): Giving a List of `u`
  61//   bez = [[0,0], [5,35], [60,-25], [80,0]];
  62//   debug_bezier(bez, N=len(bez)-1);
  63//   pts = bezier_points(bez, [0, 0.2, 0.3, 0.7, 0.8, 1]);
  64//   rainbow(pts) move($item) sphere(1.5, $fn=12);
  65// Example(2D): Giving a Range of `u`
  66//   bez = [[0,0], [5,35], [60,-25], [80,0]];
  67//   debug_bezier(bez, N=len(bez)-1);
  68//   pts = bezier_points(bez, [0:0.2:1]);
  69//   rainbow(pts) move($item) sphere(1.5, $fn=12);
  70
  71// Ugly but speed optimized code for computing bezier curves using the matrix representation
  72// See https://pomax.github.io/bezierinfo/#matrix for explanation.
  73//
  74// All of the loop unrolling makes and the use of the matrix lookup table make a big difference
  75// in the speed of execution.  For orders 10 and below this code is 10-20 times faster than
  76// the recursive code using the de Casteljau method depending on the bezier order and the
  77// number of points evaluated in one call (more points is faster).  For orders 11 and above without the
  78// lookup table or hard coded powers list the code is about twice as fast as the recursive method.
  79// Note that everything I tried to simplify or tidy this code made is slower, sometimes a lot slower.
  80function bezier_points(curve, u) =
  81    is_num(u) ? bezier_points(curve,[u])[0] :
  82    let(
  83        N = len(curve)-1,
  84        M = _bezier_matrix(N)*curve
  85    )
  86    N==0 ? [for(uval=u)[1]*M] :
  87    N==1 ? [for(uval=u)[1, uval]*M] :
  88    N==2 ? [for(uval=u)[1, uval, uval^2]*M] :
  89    N==3 ? [for(uval=u)[1, uval, uval^2, uval^3]*M] :          
  90    N==4 ? [for(uval=u)[1, uval, uval^2, uval^3, uval^4]*M] :
  91    N==5 ? [for(uval=u)[1, uval, uval^2, uval^3, uval^4, uval^5]*M] :
  92    N==6 ? [for(uval=u)[1, uval, uval^2, uval^3, uval^4, uval^5,uval^6]*M] :
  93    N==7 ? [for(uval=u)[1, uval, uval^2, uval^3, uval^4, uval^5,uval^6, uval^7]*M] :
  94    N==8 ? [for(uval=u)[1, uval, uval^2, uval^3, uval^4, uval^5,uval^6, uval^7, uval^8]*M] :
  95    N==9 ? [for(uval=u)[1, uval, uval^2, uval^3, uval^4, uval^5,uval^6, uval^7, uval^8, uval^9]*M] :
  96    N==10? [for(uval=u)[1, uval, uval^2, uval^3, uval^4, uval^5,uval^6, uval^7, uval^8, uval^9, uval^10]*M] :
  97    /* N>=11 */  [for(uval=u)[for (i=[0:1:N]) uval^i]*M];
  98
  99
 100// Not public.
 101function _signed_pascals_triangle(N,tri=[[-1]]) =
 102    len(tri)==N+1 ? tri :
 103    let(last=tri[len(tri)-1])
 104    _signed_pascals_triangle(N,concat(tri,[[-1, for(i=[0:1:len(tri)-2]) (i%2==1?-1:1)*(abs(last[i])+abs(last[i+1])),len(last)%2==0? -1:1]]));
 105
 106
 107// Not public.
 108function _compute_bezier_matrix(N) =
 109    let(tri = _signed_pascals_triangle(N))
 110    [for(i=[0:N]) concat(tri[N][i]*tri[i], repeat(0,N-i))];
 111
 112
 113// The bezier matrix, which is related to Pascal's triangle, enables nonrecursive computation
 114// of bezier points.  This method is much faster than the recursive de Casteljau method
 115// in OpenScad, but we have to precompute the matrices to reap the full benefit.
 116
 117// Not public.
 118_bezier_matrix_table = [
 119    [[1]],
 120    [[ 1, 0],
 121     [-1, 1]],
 122    [[1, 0, 0],
 123     [-2, 2, 0],
 124     [1, -2, 1]],
 125    [[ 1, 0, 0, 0],
 126     [-3, 3, 0, 0],
 127     [ 3,-6, 3, 0],
 128     [-1, 3,-3, 1]],
 129    [[ 1,  0,  0, 0, 0],
 130     [-4,  4,  0, 0, 0],
 131     [ 6,-12,  6, 0, 0],
 132     [-4, 12,-12, 4, 0],
 133     [ 1, -4,  6,-4, 1]],
 134    [[  1,  0, 0,   0, 0, 0],
 135     [ -5,  5, 0,   0, 0, 0],
 136     [ 10,-20, 10,  0, 0, 0],
 137     [-10, 30,-30, 10, 0, 0],
 138     [  5,-20, 30,-20, 5, 0],
 139     [ -1,  5,-10, 10,-5, 1]],
 140    [[  1,  0,  0,  0,  0, 0, 0],
 141     [ -6,  6,  0,  0,  0, 0, 0],
 142     [ 15,-30, 15,  0,  0, 0, 0],
 143     [-20, 60,-60, 20,  0, 0, 0],
 144     [ 15,-60, 90,-60, 15, 0, 0],
 145     [ -6, 30,-60, 60,-30, 6, 0],
 146     [  1, -6, 15,-20, 15,-6, 1]],
 147    [[  1,   0,   0,   0,  0,   0, 0, 0],
 148     [ -7,   7,   0,   0,  0,   0, 0, 0],
 149     [ 21, -42,  21,   0,  0,   0, 0, 0],
 150     [-35, 105,-105,  35,  0,   0, 0, 0],
 151     [ 35,-140, 210,-140,  35,  0, 0, 0],
 152     [-21, 105,-210, 210,-105, 21, 0, 0],
 153     [  7, -42, 105,-140, 105,-42, 7, 0],
 154     [ -1,   7, -21,  35, -35, 21,-7, 1]],
 155    [[  1,   0,   0,   0,   0,   0,  0, 0, 0],
 156     [ -8,   8,   0,   0,   0,   0,  0, 0, 0],
 157     [ 28, -56,  28,   0,   0,   0,  0, 0, 0],
 158     [-56, 168,-168,  56,   0,   0,  0, 0, 0],
 159     [ 70,-280, 420,-280,  70,   0,  0, 0, 0],
 160     [-56, 280,-560, 560,-280,  56,  0, 0, 0],
 161     [ 28,-168, 420,-560, 420,-168, 28, 0, 0],
 162     [ -8,  56,-168, 280,-280, 168,-56, 8, 0],
 163     [  1,  -8,  28, -56,  70, -56, 28,-8, 1]],
 164    [[1, 0, 0, 0, 0, 0, 0,  0, 0, 0], [-9, 9, 0, 0, 0, 0, 0, 0, 0, 0], [36, -72, 36, 0, 0, 0, 0, 0, 0, 0], [-84, 252, -252, 84, 0, 0, 0, 0, 0, 0],
 165     [126, -504, 756, -504, 126, 0, 0, 0, 0, 0], [-126, 630, -1260, 1260, -630, 126, 0, 0, 0, 0], [84, -504, 1260, -1680, 1260, -504, 84, 0, 0, 0],
 166     [-36, 252, -756, 1260, -1260, 756, -252, 36, 0, 0], [9, -72, 252, -504, 630, -504, 252, -72, 9, 0], [-1, 9, -36, 84, -126, 126, -84, 36, -9, 1]],
 167    [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [-10, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0], [45, -90, 45, 0, 0, 0, 0, 0, 0, 0, 0], [-120, 360, -360, 120, 0, 0, 0, 0, 0, 0, 0],
 168     [210, -840, 1260, -840, 210, 0, 0, 0, 0, 0, 0], [-252, 1260, -2520, 2520, -1260, 252, 0, 0, 0, 0, 0],
 169     [210, -1260, 3150, -4200, 3150, -1260, 210, 0, 0, 0, 0], [-120, 840, -2520, 4200, -4200, 2520, -840, 120, 0, 0, 0],
 170     [45, -360, 1260, -2520, 3150, -2520, 1260, -360, 45, 0, 0], [-10, 90, -360, 840, -1260, 1260, -840, 360, -90, 10, 0],
 171     [1, -10, 45, -120, 210, -252, 210, -120, 45, -10, 1]]
 172];
 173
 174
 175// Not public.
 176function _bezier_matrix(N) =
 177    N>10 ? _compute_bezier_matrix(N) :
 178    _bezier_matrix_table[N];
 179
 180
 181// Function: bezier_curve()
 182// Synopsis: Computes a specified number of points on a bezier curve.
 183// SynTags: Path
 184// Topics: Bezier Curves
 185// See Also: bezier_curve(), bezier_curvature(), bezier_tangent(), bezier_derivative(), bezier_points()
 186// Usage:
 187//   path = bezier_curve(bezier, [splinesteps], [endpoint]);
 188// Description:
 189//   Takes a list of bezier control points and generates splinesteps segments (splinesteps+1 points)
 190//   along the bezier curve they define.
 191//   Points start at the first control point and are sampled uniformly along the bezier parameter.
 192//   The endpoints of the output will be *exactly* equal to the first and last bezier control points
 193//   when endpoint is true.  If endpoint is false the sampling stops one step before the final point
 194//   of the bezier curve, but you still get the same number of (more tightly spaced) points.  
 195//   The distance between the points will *not* be equidistant.  
 196//   The degree of the bezier curve is one less than the number of points in `curve`.
 197// Arguments:
 198//   bezier = The list of control points that define the Bezier curve. 
 199//   splinesteps = The number of segments to create on the bezier curve.  Default: 16
 200//   endpoint = if false then exclude the endpoint.  Default: True
 201// Example(2D): Quadratic (Degree 2) Bezier.
 202//   bez = [[0,0], [30,30], [80,0]];
 203//   move_copies(bezier_curve(bez, 8)) sphere(r=1.5, $fn=12);
 204//   debug_bezier(bez, N=len(bez)-1);
 205// Example(2D): Cubic (Degree 3) Bezier
 206//   bez = [[0,0], [5,35], [60,-25], [80,0]];
 207//   move_copies(bezier_curve(bez, 8)) sphere(r=1.5, $fn=12);
 208//   debug_bezier(bez, N=len(bez)-1);
 209// Example(2D): Degree 4 Bezier.
 210//   bez = [[0,0], [5,15], [40,20], [60,-15], [80,0]];
 211//   move_copies(bezier_curve(bez, 8)) sphere(r=1.5, $fn=12);
 212//   debug_bezier(bez, N=len(bez)-1);
 213function bezier_curve(bezier,splinesteps=16,endpoint=true) =
 214    bezier_points(bezier, lerpn(0,1,splinesteps+1,endpoint));
 215
 216
 217// Function: bezier_derivative()
 218// Synopsis: Evaluates the derivative of the bezier curve at the given point or points.
 219// Topics: Bezier Curves
 220// See Also: bezier_curvature(), bezier_tangent(), bezier_points()
 221// Usage:
 222//   deriv = bezier_derivative(bezier, u, [order]);
 223//   derivs = bezier_derivative(bezier, LIST, [order]);
 224//   derivs = bezier_derivative(bezier, RANGE, [order]);
 225// Description:
 226//   Evaluates the derivative of the bezier curve at the given parameter value or values, `u`.  The `order` gives the order of the derivative. 
 227//   The degree of the bezier curve is one less than the number of points in `bezier`.
 228// Arguments:
 229//   bezier = The list of control points that define the Bezier curve. 
 230//   u = Parameter values for evaluating the curve, given as a single value, a list or a range.
 231//   order = The order of the derivative to return.  Default: 1 (for the first derivative)
 232function bezier_derivative(bezier, u, order=1) =
 233    assert(is_int(order) && order>=0)
 234    order==0? bezier_points(bezier, u) : let(
 235        N = len(bezier) - 1,
 236        dpts = N * deltas(bezier)
 237    ) order==1? bezier_points(dpts, u) :
 238    bezier_derivative(dpts, u, order-1);
 239
 240
 241
 242// Function: bezier_tangent()
 243// Synopsis: Calculates unit tangent vectors along the bezier curve at one or more given positions.
 244// Topics: Bezier Curves
 245// See Also: bezier_curvature(), bezier_derivative(), bezier_points()
 246// Usage:
 247//   tanvec = bezier_tangent(bezier, u);
 248//   tanvecs = bezier_tangent(bezier, LIST);
 249//   tanvecs = bezier_tangent(bezier, RANGE);
 250// Description:
 251//   Returns the unit tangent vector at the given parameter values on a bezier curve with control points `bezier`.
 252// Arguments:
 253//   bezier = The list of control points that define the Bezier curve. 
 254//   u = Parameter values for evaluating the curve, given as a single value, a list or a range.
 255function bezier_tangent(bezier, u) =
 256    let(
 257        res = bezier_derivative(bezier, u)
 258    ) is_vector(res)? unit(res) :
 259    [for (v=res) unit(v)];
 260
 261
 262
 263// Function: bezier_curvature()
 264// Synopsis: Returns the curvature at one or more given positions along a bezier curve.
 265// Topics: Bezier Curves
 266// See Also: bezier_tangent(), bezier_derivative(), bezier_points()
 267// Usage:
 268//   crv = bezier_curvature(curve, u);
 269//   crvlist = bezier_curvature(curve, LIST);
 270//   crvlist = bezier_curvature(curve, RANGE);
 271// Description:
 272//   Returns the curvature value for the given parameters `u` on the bezier curve with control points `bezier`. 
 273//   The curvature is the inverse of the radius of the tangent circle at the given point.
 274//   Thus, the tighter the curve, the larger the curvature value.  Curvature will be 0 for
 275//   a position with no curvature, since 1/0 is not a number.
 276// Arguments:
 277//   bezier = The list of control points that define the Bezier curve.
 278//   u = Parameter values for evaluating the curve, given as a single value, a list or a range.
 279function bezier_curvature(bezier, u) =
 280    is_num(u) ? bezier_curvature(bezier,[u])[0] :
 281    let(
 282        d1 = bezier_derivative(bezier, u, 1),
 283        d2 = bezier_derivative(bezier, u, 2)
 284    ) [
 285        for(i=idx(d1))
 286        sqrt(
 287            sqr(norm(d1[i])*norm(d2[i])) -
 288            sqr(d1[i]*d2[i])
 289        ) / pow(norm(d1[i]),3)
 290    ];
 291
 292
 293
 294// Function: bezier_closest_point()
 295// Synopsis: Finds the closest position on a bezier curve to a given point.
 296// Topics: Bezier Curves
 297// See Also: bezier_points()
 298// Usage:
 299//   u = bezier_closest_point(bezier, pt, [max_err]);
 300// Description:
 301//   Finds the closest part of the given bezier curve to point `pt`.
 302//   The degree of the curve, N, is one less than the number of points in `curve`.
 303//   Returns `u` for the closest position on the bezier curve to the given point `pt`.
 304// Arguments:
 305//   bezier = The list of control points that define the Bezier curve. 
 306//   pt = The point to find the closest curve point to.
 307//   max_err = The maximum allowed error when approximating the closest approach.
 308// Example(2D):
 309//   pt = [40,15];
 310//   bez = [[0,0], [20,40], [60,-25], [80,0]];
 311//   u = bezier_closest_point(bez, pt);
 312//   debug_bezier(bez, N=len(bez)-1);
 313//   color("red") translate(pt) sphere(r=1);
 314//   color("blue") translate(bezier_points(bez,u)) sphere(r=1);
 315function bezier_closest_point(bezier, pt, max_err=0.01, u=0, end_u=1) =
 316    let(
 317        steps = len(bezier)*3,
 318        uvals = [u, for (i=[0:1:steps]) (end_u-u)*(i/steps)+u, end_u],
 319        path = bezier_points(bezier,uvals),
 320        minima_ranges = [
 321            for (i = [1:1:len(uvals)-2]) let(
 322                d1 = norm(path[i-1]-pt),
 323                d2 = norm(path[i  ]-pt),
 324                d3 = norm(path[i+1]-pt)
 325            ) if (d2<=d1 && d2<=d3) [uvals[i-1],uvals[i+1]]
 326        ]
 327    ) len(minima_ranges)>1? (
 328        let(
 329            min_us = [
 330                for (minima = minima_ranges)
 331                    bezier_closest_point(bezier, pt, max_err=max_err, u=minima.x, end_u=minima.y)
 332            ],
 333            dists = [for (v=min_us) norm(bezier_points(bezier,v)-pt)],
 334            min_i = min_index(dists)
 335        ) min_us[min_i]
 336    ) : let(
 337        minima = minima_ranges[0],
 338        pp = bezier_points(bezier, minima),
 339        err = norm(pp[1]-pp[0])
 340    ) err<max_err? mean(minima) :
 341    bezier_closest_point(bezier, pt, max_err=max_err, u=minima[0], end_u=minima[1]);
 342
 343
 344// Function: bezier_length()
 345// Synopsis: Approximate the length of part of a bezier curve.
 346// Topics: Bezier Curves
 347// See Also: bezier_points()
 348// Usage:
 349//   pathlen = bezier_length(bezier, [start_u], [end_u], [max_deflect]);
 350// Description:
 351//   Approximates the length of the portion of the bezier curve between start_u and end_u.
 352// Arguments:
 353//   bezier = The list of control points that define the Bezier curve. 
 354//   start_u = The Bezier parameter to start measuring measuring from.  Between 0 and 1.
 355//   end_u = The Bezier parameter to end measuring at.  Between 0 and 1.  Greater than start_u.
 356//   max_deflect = The largest amount of deflection from the true curve to allow for approximation.
 357// Example:
 358//   bez = [[0,0], [5,35], [60,-25], [80,0]];
 359//   echo(bezier_length(bez));
 360function bezier_length(bezier, start_u=0, end_u=1, max_deflect=0.01) =
 361    let(
 362        segs = len(bezier) * 2,
 363        uvals = lerpn(start_u, end_u, segs+1),
 364        path = bezier_points(bezier,uvals),
 365        defl = max([
 366            for (i=idx(path,e=-3)) let(
 367                mp = (path[i] + path[i+2]) / 2
 368            ) norm(path[i+1] - mp)
 369        ]),
 370        mid_u = lerp(start_u, end_u, 0.5)
 371    )
 372    defl <= max_deflect? path_length(path) :
 373    sum([
 374        for (i=[0:1:segs-1]) let(
 375            su = lerp(start_u, end_u, i/segs),
 376            eu = lerp(start_u, end_u, (i+1)/segs)
 377        ) bezier_length(bezier, su, eu, max_deflect)
 378    ]);
 379
 380
 381
 382// Function: bezier_line_intersection()
 383// Synopsis: Calculates where a bezier curve intersects a line.
 384// Topics: Bezier Curves, Geometry, Intersection
 385// See Also: bezier_points(), bezier_length(), bezier_closest_point()
 386// Usage: 
 387//   u = bezier_line_intersection(bezier, line);
 388// Description:
 389//   Finds the intersection points of the 2D Bezier curve with control points `bezier` and the given line, specified as a pair of points.  
 390//   Returns the intersection as a list of `u` values for the Bezier.  
 391// Arguments:
 392//   bezier = The list of control points that define a 2D Bezier curve. 
 393//   line = a list of two distinct 2d points defining a line
 394function bezier_line_intersection(bezier, line) =
 395    assert(is_path(bezier,2), "The input ´bezier´ must be a 2d bezier")
 396    assert(_valid_line(line,2), "The input `line` is not a valid 2d line")
 397    let( 
 398        a = _bezier_matrix(len(bezier)-1)*bezier, // bezier algebraic coeffs. 
 399        n = [-line[1].y+line[0].y, line[1].x-line[0].x], // line normal
 400        q = [for(i=[len(a)-1:-1:1]) a[i]*n, (a[0]-line[0])*n] // bezier SDF to line
 401    )
 402    [for(u=real_roots(q)) if (u>=0 && u<=1) u];
 403
 404
 405
 406
 407// Section: Bezier Path Functions
 408//   To contruct more complicated curves you can connect a sequence of Bezier curves end to end.  
 409//   A Bezier path is a flattened list of control points that, along with the degree, represents such a sequence of bezier curves where all of the curves have the same degree.
 410//   A Bezier path looks like a regular path, since it is just a list of points, but it is not a regular path.  Use {{bezpath_curve()}} to convert a Bezier path to a regular path.
 411//   We interpret a degree N Bezier path as groups of N+1 control points that
 412//   share endpoints, so they overlap by one point.  So if you have an order 3 bezier path `[p0,p1,p2,p3,p4,p5,p6]` then the first
 413//   Bezier curve control point set is `[p0,p1,p2,p3]` and the second one is `[p3,p4,p5,p6]`.  The endpoint, `p3`, is shared between the control point sets.
 414//   The Bezier degree, which must be known to interpret the Bezier path, defaults to 3. 
 415
 416
 417// Function: bezpath_points()
 418// Synopsis: Computes one or more specified points along a bezier path.
 419// SynTags: Path
 420// Topics: Bezier Paths
 421// See Also: bezier_points(), bezier_curve()
 422// Usage:
 423//   pt = bezpath_points(bezpath, curveind, u, [N]);
 424//   ptlist = bezpath_points(bezpath, curveind, LIST, [N]);
 425//   path = bezpath_points(bezpath, curveind, RANGE, [N]);
 426// Description:
 427//   Extracts from the Bezier path `bezpath` the control points for the Bezier curve whose index is `curveind` and
 428//   computes the point or points on the corresponding Bezier curve specified by `u`.  If `curveind` is zero you
 429//   get the first curve.  The number of curves is `(len(bezpath)-1)/N` so the maximum index is that number minus one.  
 430// Arguments:
 431//   bezpath = A Bezier path path to approximate.
 432//   curveind = Curve number along the path.  
 433//   u = Parameter values for evaluating the curve, given as a single value, a list or a range.
 434//   N = The degree of the Bezier path curves.  Default: 3
 435function bezpath_points(bezpath, curveind, u, N=3) =
 436    bezier_points(select(bezpath,curveind*N,(curveind+1)*N), u);
 437
 438
 439// Function: bezpath_curve()
 440// Synopsis: Converts bezier path into a path of points. 
 441// SynTags: Path
 442// Topics: Bezier Paths
 443// See Also: bezier_points(), bezier_curve(), bezpath_points()
 444// Usage:
 445//   path = bezpath_curve(bezpath, [splinesteps], [N], [endpoint])
 446// Description:
 447//   Computes a number of uniformly distributed points along a bezier path.
 448// Arguments:
 449//   bezpath = A bezier path to approximate.
 450//   splinesteps = Number of straight lines to split each bezier curve into. default=16
 451//   N = The degree of the bezier curves.  Cubic beziers have N=3.  Default: 3
 452//   endpoint = If true, include the very last point of the bezier path.  Default: true
 453// Example(2D):
 454//   bez = [
 455//       [0,0], [-5,30],
 456//       [20,60], [50,50], [110,30],
 457//       [60,25], [70,0], [80,-25],
 458//       [80,-50], [50,-50]
 459//   ];
 460//   path = bezpath_curve(bez);
 461//   stroke(path,dots=true,dots_color="red");
 462function bezpath_curve(bezpath, splinesteps=16, N=3, endpoint=true) =
 463    assert(is_path(bezpath))
 464    assert(is_int(N))
 465    assert(is_int(splinesteps) && splinesteps>0)
 466    assert(len(bezpath)%N == 1, str("A degree ",N," bezier path should have a multiple of ",N," points in it, plus 1."))
 467    let(
 468        segs = (len(bezpath)-1) / N,
 469        step = 1 / splinesteps
 470    ) [
 471        for (seg = [0:1:segs-1])
 472            each bezier_points(select(bezpath, seg*N, (seg+1)*N), [0:step:1-step/2]),
 473        if (endpoint) last(bezpath)
 474    ];
 475
 476
 477// Function: bezpath_closest_point()
 478// Synopsis: Finds the closest point on a bezier path to a given point.
 479// Topics: Bezier Paths
 480// See Also: bezpath_points(), bezpath_curve(), bezier_points(), bezier_curve(), bezier_closest_point()
 481// Usage:
 482//   res = bezpath_closest_point(bezpath, pt, [N], [max_err]);
 483// Description:
 484//   Finds an approximation to the closest part of the given bezier path to point `pt`.
 485//   Returns [segnum, u] for the closest position on the bezier path to the given point `pt`.
 486// Arguments:
 487//   bezpath = A bezier path to approximate.
 488//   pt = The point to find the closest curve point to.
 489//   N = The degree of the bezier curves.  Cubic beziers have N=3.  Default: 3
 490//   max_err = The maximum allowed error when approximating the closest approach.
 491// Example(2D):
 492//   pt = [100,0];
 493//   bez = [[0,0], [20,40], [60,-25], [80,0],
 494//          [100,25], [140,25], [160,0]];
 495//   pos = bezpath_closest_point(bez, pt);
 496//   xy = bezpath_points(bez,pos[0],pos[1]);
 497//   debug_bezier(bez, N=3);
 498//   color("red") translate(pt) sphere(r=1);
 499//   color("blue") translate(xy) sphere(r=1);
 500function bezpath_closest_point(bezpath, pt, N=3, max_err=0.01, seg=0, min_seg=undef, min_u=undef, min_dist=undef) =
 501    assert(is_vector(pt))
 502    assert(is_int(N))
 503    assert(is_num(max_err))
 504    assert(len(bezpath)%N == 1, str("A degree ",N," bezier path shound have a multiple of ",N," points in it, plus 1."))
 505    let(curve = select(bezpath,seg*N,(seg+1)*N))
 506    (seg*N+1 >= len(bezpath))? (
 507        let(curve = select(bezpath, min_seg*N, (min_seg+1)*N))
 508        [min_seg, bezier_closest_point(curve, pt, max_err=max_err)]
 509    ) : (
 510        let(
 511            curve = select(bezpath,seg*N,(seg+1)*N),
 512            u = bezier_closest_point(curve, pt, max_err=0.05),
 513            dist = norm(bezier_points(curve, u)-pt),
 514            mseg = (min_dist==undef || dist<min_dist)? seg : min_seg,
 515            mdist = (min_dist==undef || dist<min_dist)? dist : min_dist,
 516            mu = (min_dist==undef || dist<min_dist)? u : min_u
 517        )
 518        bezpath_closest_point(bezpath, pt, N, max_err, seg+1, mseg, mu, mdist)
 519    );
 520
 521
 522
 523// Function: bezpath_length()
 524// Synopsis: Approximate the length of a bezier path.
 525// Topics: Bezier Paths
 526// See Also: bezier_points(), bezier_curve(), bezier_length()
 527// Usage:
 528//   plen = bezpath_length(path, [N], [max_deflect]);
 529// Description:
 530//   Approximates the length of the bezier path.
 531// Arguments:
 532//   path = A bezier path to approximate.
 533//   N = The degree of the bezier curves.  Cubic beziers have N=3.  Default: 3
 534//   max_deflect = The largest amount of deflection from the true curve to allow for approximation.
 535function bezpath_length(bezpath, N=3, max_deflect=0.001) =
 536    assert(is_int(N))
 537    assert(is_num(max_deflect))
 538    assert(len(bezpath)%N == 1, str("A degree ",N," bezier path shound have a multiple of ",N," points in it, plus 1."))
 539    sum([
 540        for (seg=[0:1:(len(bezpath)-1)/N-1]) (
 541            bezier_length(
 542                select(bezpath, seg*N, (seg+1)*N),
 543                max_deflect=max_deflect
 544            )
 545        )
 546    ]);
 547
 548
 549
 550// Function: path_to_bezpath()
 551// Synopsis: Generates a bezier path that passes through all points in a given linear path.
 552// SynTags: Path
 553// Topics: Bezier Paths, Rounding
 554// See Also: path_tangents()
 555// Usage:
 556//   bezpath = path_to_bezpath(path, [closed], [tangents], [uniform], [size=]|[relsize=]);
 557// Description:
 558//   Given a 2d or 3d input path and optional list of tangent vectors, computes a cubic (degree 3) bezier
 559//   path that passes through every point on the input path and matches the tangent vectors.  If you do
 560//   not supply the tangent it will be computed using `path_tangents()`.  If the path is closed specify this
 561//   by setting `closed=true`.  The size or relsize parameter determines how far the curve can deviate from
 562//   the input path.  In the case where the curve has a single hump, the size specifies the exact distance
 563//   between the specified path and the bezier.  If you give relsize then it is relative to the segment
 564//   length (e.g. 0.05 means 5% of the segment length).  In 2d when the bezier curve makes an S-curve
 565//   the size parameter specifies the sum of the deviations of the two peaks of the curve.  In 3-space
 566//   the bezier curve may have three extrema: two maxima and one minimum.  In this case the size specifies
 567//   the sum of the maxima minus the minimum.  If you do not supply the tangents then they are computed
 568//   using `path_tangents()` with `uniform=false` by default.  Tangents computed on non-uniform data tend
 569//   to display overshoots.  See `smooth_path()` for examples.
 570// Arguments:
 571//   path = 2D or 3D point list or 1-region that the curve must pass through
 572//   closed = true if the curve is closed .  Default: false
 573//   tangents = tangents constraining curve direction at each point
 574//   uniform = set to true to compute tangents with uniform=true.  Default: false
 575//   ---
 576//   size = absolute size specification for the curve, a number or vector
 577//   relsize = relative size specification for the curve, a number or vector.  Default: 0.1. 
 578function path_to_bezpath(path, closed, tangents, uniform=false, size, relsize) =
 579    is_1region(path) ? path_to_bezpath(path[0], default(closed,true), tangents, uniform, size, relsize) :
 580    let(closed=default(closed,false))
 581    assert(is_bool(closed))
 582    assert(is_bool(uniform))
 583    assert(num_defined([size,relsize])<=1, "Can't define both size and relsize")
 584    assert(is_path(path,[2,3]),"Input path is not a valid 2d or 3d path")
 585    assert(is_undef(tangents) || is_path(tangents,[2,3]),"Tangents must be a 2d or 3d path")
 586    assert(is_undef(tangents) || len(path)==len(tangents), "Input tangents must be the same length as the input path")
 587    let(
 588        curvesize = first_defined([size,relsize,0.1]),
 589        relative = is_undef(size),
 590        lastpt = len(path) - (closed?0:1)
 591    )
 592    assert(is_num(curvesize) || len(curvesize)==lastpt, str("Size or relsize must have length ",lastpt))
 593    let(
 594        sizevect = is_num(curvesize) ? repeat(curvesize, lastpt) : curvesize,
 595        tangents = is_def(tangents) ? [for(t=tangents) let(n=norm(t)) assert(!approx(n,0),"Zero tangent vector") t/n] :
 596                                      path_tangents(path, uniform=uniform, closed=closed)
 597    )
 598    assert(min(sizevect)>0, "Size and relsize must be greater than zero")
 599    [
 600        for(i=[0:1:lastpt-1])
 601            let(
 602                first = path[i],
 603                second = select(path,i+1),
 604                seglength = norm(second-first),
 605                dummy = assert(seglength>0, str("Path segment has zero length from index ",i," to ",i+1)),
 606                segdir = (second-first)/seglength,
 607                tangent1 = tangents[i],
 608                tangent2 = -select(tangents,i+1),                        // Need this to point backwards, in direction of the curve
 609                parallel = abs(tangent1*segdir) + abs(tangent2*segdir), // Total component of tangents parallel to the segment
 610                Lmax = seglength/parallel,    // May be infinity
 611                size = relative ? sizevect[i]*seglength : sizevect[i],
 612                normal1 = tangent1-(tangent1*segdir)*segdir,   // Components of the tangents orthogonal to the segment
 613                normal2 = tangent2-(tangent2*segdir)*segdir,
 614                p = [ [-3 ,6,-3 ],                   // polynomial in power form
 615                      [ 7,-9, 2 ],
 616                      [-5, 3, 0 ],
 617                      [ 1, 0, 0 ] ]*[normal1*normal1, normal1*normal2, normal2*normal2],
 618                uextreme = approx(norm(p),0) ? []
 619                                             : [for(root = real_roots(p)) if (root>0 && root<1) root],
 620                distlist = [for(d=bezier_points([normal1*0, normal1, normal2, normal2*0], uextreme)) norm(d)],
 621                scale = len(distlist)==0 ? 0 :
 622                        len(distlist)==1 ? distlist[0]
 623                                         : sum(distlist) - 2*min(distlist),
 624                Ldesired = size/scale,   // This will be infinity when the polynomial is zero
 625                L = min(Lmax, Ldesired)
 626            )
 627            each [
 628                  first, 
 629                  first + L*tangent1,
 630                  second + L*tangent2 
 631                 ],
 632        select(path,lastpt)
 633    ];
 634
 635
 636
 637
 638// Function: bezpath_close_to_axis()
 639// Synopsis: Closes a 2D bezier path to the specified axis.
 640// SynTags: Path
 641// Topics: Bezier Paths
 642// See Also: bezpath_offset()
 643// Usage:
 644//   bezpath = bezpath_close_to_axis(bezpath, [axis], [N]);
 645// Description:
 646//   Takes a 2D bezier path and closes it to the specified axis.
 647// Arguments:
 648//   bezpath = The 2D bezier path to close to the axis.
 649//   axis = The axis to close to, "X", or "Y".  Default: "X"
 650//   N = The degree of the bezier curves.  Cubic beziers have N=3.  Default: 3
 651// Example(2D):
 652//   bez = [[50,30], [40,10], [10,50], [0,30],
 653//          [-10, 10], [-30,10], [-50,20]];
 654//   closed = bezpath_close_to_axis(bez);
 655//   debug_bezier(closed);
 656// Example(2D):
 657//   bez = [[30,50], [10,40], [50,10], [30,0],
 658//          [10, -10], [10,-30], [20,-50]];
 659//   closed = bezpath_close_to_axis(bez, axis="Y");
 660//   debug_bezier(closed);
 661function bezpath_close_to_axis(bezpath, axis="X", N=3) =
 662    assert(is_path(bezpath,2), "bezpath_close_to_axis() can only work on 2D bezier paths.")
 663    assert(is_int(N))
 664    assert(len(bezpath)%N == 1, str("A degree ",N," bezier path shound have a multiple of ",N," points in it, plus 1."))
 665    let(
 666        sp = bezpath[0],
 667        ep = last(bezpath)
 668    ) (axis=="X")? concat(
 669        lerpn([sp.x,0], sp, N, false),
 670        list_head(bezpath),
 671        lerpn(ep, [ep.x,0], N, false),
 672        lerpn([ep.x,0], [sp.x,0], N+1)
 673    ) : (axis=="Y")? concat(
 674        lerpn([0,sp.y], sp, N, false),
 675        list_head(bezpath),
 676        lerpn(ep, [0,ep.y], N, false),
 677        lerpn([0,ep.y], [0,sp.y], N+1)
 678    ) : (
 679        assert(in_list(axis, ["X","Y"]))
 680    );
 681
 682
 683// Function: bezpath_offset()
 684// Synopsis: Forms a closed bezier path loop with a translated and reversed copy of itself.
 685// SynTags: Path
 686// Topics: Bezier Paths
 687// See Also: bezpath_close_to_axis()
 688// Usage:
 689//   bezpath = bezpath_offset(offset, bezier, [N]);
 690// Description:
 691//   Takes a 2D bezier path and closes it with a matching reversed path that is offset by the given `offset` [X,Y] distance.
 692// Arguments:
 693//   offset = Amount to offset second path by.
 694//   bezier = The 2D bezier path.
 695//   N = The degree of the bezier curves.  Cubic beziers have N=3.  Default: 3
 696// Example(2D):
 697//   bez = [[50,30], [40,10], [10,50], [0,30], [-10, 10], [-30,10], [-50,20]];
 698//   closed = bezpath_offset([0,-5], bez);
 699//   debug_bezier(closed);
 700// Example(2D):
 701//   bez = [[30,50], [10,40], [50,10], [30,0], [10, -10], [10,-30], [20,-50]];
 702//   closed = bezpath_offset([-5,0], bez);
 703//   debug_bezier(closed);
 704function bezpath_offset(offset, bezier, N=3) =
 705    assert(is_vector(offset,2))
 706    assert(is_path(bezier,2), "bezpath_offset() can only work on 2D bezier paths.")
 707    assert(is_int(N))
 708    assert(len(bezier)%N == 1, str("A degree ",N," bezier path shound have a multiple of ",N," points in it, plus 1."))
 709    let(
 710        backbez = reverse([ for (pt = bezier) pt+offset ]),
 711        bezend = len(bezier)-1
 712    ) concat(
 713        list_head(bezier),
 714        lerpn(bezier[bezend], backbez[0], N, false),
 715        list_head(backbez),
 716        lerpn(backbez[bezend], bezier[0], N+1)
 717    );
 718
 719
 720
 721// Section: Cubic Bezier Path Construction
 722
 723// Function: bez_begin()
 724// Synopsis: Calculates starting bezier path control points.
 725// Topics: Bezier Paths
 726// See Also: bez_tang(), bez_joint(), bez_end()
 727// Usage:
 728//   pts = bez_begin(pt, a, r, [p=]);
 729//   pts = bez_begin(pt, VECTOR, [r], [p=]);
 730// Description:
 731//   This is used to create the first endpoint and control point of a cubic bezier path.
 732// Arguments:
 733//   pt = The starting endpoint for the bezier path.
 734//   a = If given a scalar, specifies the theta (XY plane) angle in degrees from X+.  If given a vector, specifies the direction and possibly distance of the first control point.
 735//   r = Specifies the distance of the control point from the endpoint `pt`.
 736//   ---
 737//   p = If given, specifies the number of degrees away from the Z+ axis.
 738// Example(2D): 2D Bezier Path by Angle
 739//   bezpath = flatten([
 740//       bez_begin([-50,  0],  45,20),
 741//       bez_tang ([  0,  0],-135,20),
 742//       bez_joint([ 20,-25], 135, 90, 10, 15),
 743//       bez_end  ([ 50,  0], -90,20),
 744//   ]);
 745//   debug_bezier(bezpath);
 746// Example(2D): 2D Bezier Path by Vector
 747//   bezpath = flatten([
 748//       bez_begin([-50,0],[0,-20]),
 749//       bez_tang ([-10,0],[0,-20]),
 750//       bez_joint([ 20,-25], [-10,10], [0,15]),
 751//       bez_end  ([ 50,0],[0, 20]),
 752//   ]);
 753//   debug_bezier(bezpath);
 754// Example(2D): 2D Bezier Path by Vector and Distance
 755//   bezpath = flatten([
 756//       bez_begin([-30,0],FWD, 30),
 757//       bez_tang ([  0,0],FWD, 30),
 758//       bez_joint([ 20,-25], 135, 90, 10, 15),
 759//       bez_end  ([ 30,0],BACK,30),
 760//   ]);
 761//   debug_bezier(bezpath);
 762// Example(3D,FlatSpin,VPD=200): 3D Bezier Path by Angle
 763//   bezpath = flatten([
 764//       bez_begin([-30,0,0],90,20,p=135),
 765//       bez_tang ([  0,0,0],-90,20,p=135),
 766//       bez_joint([20,-25,0], 135, 90, 15, 10, p1=135, p2=45),
 767//       bez_end  ([ 30,0,0],-90,20,p=45),
 768//   ]);
 769//   debug_bezier(bezpath);
 770// Example(3D,FlatSpin,VPD=225): 3D Bezier Path by Vector
 771//   bezpath = flatten([
 772//       bez_begin([-30,0,0],[0,-20, 20]),
 773//       bez_tang ([  0,0,0],[0,-20,-20]),
 774//       bez_joint([20,-25,0],[0,10,-10],[0,15,15]),
 775//       bez_end  ([ 30,0,0],[0,-20,-20]),
 776//   ]);
 777//   debug_bezier(bezpath);
 778// Example(3D,FlatSpin,VPD=225): 3D Bezier Path by Vector and Distance
 779//   bezpath = flatten([
 780//       bez_begin([-30,0,0],FWD, 20),
 781//       bez_tang ([  0,0,0],DOWN,20),
 782//       bez_joint([20,-25,0],LEFT,DOWN,r1=20,r2=15),
 783//       bez_end  ([ 30,0,0],DOWN,20),
 784//   ]);
 785//   debug_bezier(bezpath);
 786function bez_begin(pt,a,r,p) =
 787    assert(is_finite(r) || is_vector(a))
 788    assert(len(pt)==3 || is_undef(p))
 789    is_vector(a)? [pt, pt+(is_undef(r)? a : r*unit(a))] :
 790    is_finite(a)? [pt, pt+spherical_to_xyz(r,a,default(p,90))] :
 791    assert(false, "Bad arguments.");
 792
 793
 794// Function: bez_tang()
 795// Synopsis: Calculates control points for a smooth bezier path joint.
 796// Topics: Bezier Paths
 797// See Also: bez_begin(), bez_joint(), bez_end()
 798// Usage:
 799//   pts = bez_tang(pt, a, r1, r2, [p=]);
 800//   pts = bez_tang(pt, VECTOR, [r1], [r2], [p=]);
 801// Description:
 802//   This creates a smooth joint in a cubic bezier path.  It creates three points, being the
 803//   approaching control point, the fixed bezier control point, and the departing control
 804//   point.  The two control points will be collinear with the fixed point, making for a
 805//   smooth bezier curve at the fixed point. See {{bez_begin()}} for examples.
 806// Arguments:
 807//   pt = The fixed point for the bezier path.
 808//   a = If given a scalar, specifies the theta (XY plane) angle in degrees from X+.  If given a vector, specifies the direction and possibly distance of the departing control point.
 809//   r1 = Specifies the distance of the approching control point from the fixed point.  Overrides the distance component of the vector if `a` contains a vector.
 810//   r2 = Specifies the distance of the departing control point from the fixed point.  Overrides the distance component of the vector if `a` contains a vector.  If `r1` is given and `r2` is not, uses the value of `r1` for `r2`.
 811//   ---
 812//   p = If given, specifies the number of degrees away from the Z+ axis.
 813function bez_tang(pt,a,r1,r2,p) =
 814    assert(is_finite(r1) || is_vector(a))
 815    assert(len(pt)==3 || is_undef(p))
 816    let(
 817        r1 = is_num(r1)? r1 : norm(a),
 818        r2 = default(r2,r1),
 819        p = default(p, 90)
 820    )
 821    is_vector(a)? [pt-r1*unit(a), pt, pt+r2*unit(a)] :
 822    is_finite(a)? [
 823        pt-spherical_to_xyz(r1,a,p),
 824        pt,
 825        pt+spherical_to_xyz(r2,a,p)
 826    ] :
 827    assert(false, "Bad arguments.");
 828
 829
 830// Function: bez_joint()
 831// Synopsis: Calculates control points for a disjointed corner bezier path joint.
 832// Topics: Bezier Paths
 833// See Also: bez_begin(), bez_tang(), bez_end()
 834// Usage:
 835//   pts = bez_joint(pt, a1, a2, r1, r2, [p1=], [p2=]);
 836//   pts = bez_joint(pt, VEC1, VEC2, [r1=], [r2=], [p1=], [p2=]);
 837// Description:
 838//   This creates a disjoint corner joint in a cubic bezier path.  It creates three points, being
 839//   the aproaching control point, the fixed bezier control point, and the departing control point.
 840//   The two control points can be directed in different arbitrary directions from the fixed bezier
 841//   point. See {{bez_begin()}} for examples.
 842// Arguments:
 843//   pt = The fixed point for the bezier path.
 844//   a1 = If given a scalar, specifies the theta (XY plane) angle in degrees from X+.  If given a vector, specifies the direction and possibly distance of the approaching control point.
 845//   a2 = If given a scalar, specifies the theta (XY plane) angle in degrees from X+.  If given a vector, specifies the direction and possibly distance of the departing control point.
 846//   r1 = Specifies the distance of the approching control point from the fixed point.  Overrides the distance component of the vector if `a1` contains a vector.
 847//   r2 = Specifies the distance of the departing control point from the fixed point.  Overrides the distance component of the vector if `a2` contains a vector.
 848//   ---
 849//   p1 = If given, specifies the number of degrees away from the Z+ axis of the approaching control point.
 850//   p2 = If given, specifies the number of degrees away from the Z+ axis of the departing control point.
 851function bez_joint(pt,a1,a2,r1,r2,p1,p2) =
 852    assert(is_finite(r1) || is_vector(a1))
 853    assert(is_finite(r2) || is_vector(a2))
 854    assert(len(pt)==3 || (is_undef(p1) && is_undef(p2)))
 855    let(
 856        r1 = is_num(r1)? r1 : norm(a1),
 857        r2 = is_num(r2)? r2 : norm(a2),
 858        p1 = default(p1, 90),
 859        p2 = default(p2, 90)
 860    ) [
 861        if (is_vector(a1)) (pt+r1*unit(a1))
 862        else if (is_finite(a1)) (pt+spherical_to_xyz(r1,a1,p1))
 863        else assert(false, "Bad Arguments"),
 864        pt,
 865        if (is_vector(a2)) (pt+r2*unit(a2))
 866        else if (is_finite(a2)) (pt+spherical_to_xyz(r2,a2,p2))
 867        else assert(false, "Bad Arguments")
 868    ];
 869
 870
 871// Function: bez_end()
 872// Synopsis: Calculates ending bezier path control points.
 873// Topics: Bezier Paths
 874// See Also: bez_tang(), bez_joint(), bez_end()
 875// Usage:
 876//   pts = bez_end(pt, a, r, [p=]);
 877//   pts = bez_end(pt, VECTOR, [r], [p=]);
 878// Description:
 879//   This is used to create the approaching control point, and the endpoint of a cubic bezier path.
 880//   See {{bez_begin()}} for examples.
 881// Arguments:
 882//   pt = The starting endpoint for the bezier path.
 883//   a = If given a scalar, specifies the theta (XY plane) angle in degrees from X+.  If given a vector, specifies the direction and possibly distance of the first control point.
 884//   r = Specifies the distance of the control point from the endpoint `pt`.
 885//   p = If given, specifies the number of degrees away from the Z+ axis.
 886function bez_end(pt,a,r,p) =
 887    assert(is_finite(r) || is_vector(a))
 888    assert(len(pt)==3 || is_undef(p))
 889    is_vector(a)? [pt+(is_undef(r)? a : r*unit(a)), pt] :
 890    is_finite(a)? [pt+spherical_to_xyz(r,a,default(p,90)), pt] :
 891    assert(false, "Bad arguments.");
 892
 893
 894
 895// Section: Bezier Surfaces
 896
 897
 898// Function: is_bezier_patch()
 899// Synopsis: Returns true if the given item is a bezier patch.
 900// Topics: Bezier Patches, Type Checking
 901// Usage:
 902//   bool = is_bezier_patch(x);
 903// Description:
 904//   Returns true if the given item is a bezier patch. (a 2D array of 3D points.)
 905// Arguments:
 906//   x = The value to check the type of.
 907function is_bezier_patch(x) =
 908    is_list(x) && is_list(x[0]) && is_vector(x[0][0]) && len(x[0]) == len(x[len(x)-1]);  
 909
 910
 911// Function: bezier_patch_flat()
 912// Synopsis: Creates a flat bezier patch.
 913// Topics: Bezier Patches
 914// See Also: bezier_patch_points()
 915// Usage:
 916//   patch = bezier_patch_flat(size, [N=], [spin=], [orient=], [trans=]);
 917// Description:
 918//   Returns a flat rectangular bezier patch of degree `N`, centered on the XY plane.
 919// Arguments:
 920//   size = scalar or 2-vector giving the X and Y dimensions of the patch. 
 921//   ---
 922//   N = Degree of the patch to generate.  Since this is flat, a degree of 1 should usually be sufficient.  Default: 1
 923//   orient = A direction vector.  Point the patch normal in this direction.  
 924//   spin = Spin angle to apply to the patch
 925//   trans = Amount to translate patch, after orient and spin. 
 926// Example(3D):
 927//   patch = bezier_patch_flat(size=[100,100]);
 928//   debug_bezier_patches([patch], size=1, showcps=true);
 929function bezier_patch_flat(size, N=1, spin=0, orient=UP, trans=[0,0,0]) =
 930    assert(N>0)
 931    let(size = force_list(size,2))
 932    assert(is_vector(size,2))
 933    let(
 934        patch = [
 935            for (x=[0:1:N]) [
 936                for (y=[0:1:N])
 937                v_mul(point3d(size), [x/N-0.5, 0.5-y/N, 0])
 938            ]
 939        ],
 940        m = move(trans) * rot(a=spin, from=UP, to=orient)
 941    ) [for (row=patch) apply(m, row)];
 942
 943
 944
 945// Function: bezier_patch_reverse()
 946// Synopsis: Reverses the orientation of a bezier patch.
 947// Topics: Bezier Patches
 948// See Also: bezier_patch_points(), bezier_patch_flat()
 949// Usage:
 950//   rpatch = bezier_patch_reverse(patch);
 951// Description:
 952//   Reverses the patch, so that the faces generated from it are flipped back to front.
 953// Arguments:
 954//   patch = The patch to reverse.
 955function bezier_patch_reverse(patch) =
 956    [for (row=patch) reverse(row)];
 957
 958
 959// Function: bezier_patch_points()
 960// Synopsis: Computes one or more specified points across a bezier surface patch.
 961// Topics: Bezier Patches
 962// See Also: bezier_patch_normals(), bezier_points(), bezier_curve(), bezpath_curve()
 963// Usage:
 964//   pt = bezier_patch_points(patch, u, v);
 965//   ptgrid = bezier_patch_points(patch, LIST, LIST);
 966//   ptgrid = bezier_patch_points(patch, RANGE, RANGE);
 967// Description:
 968//   Sample a bezier patch on a listed point set.  The bezier patch must be a rectangular array of
 969//   points, and it will be sampled at all the (u,v) pairs that you specify.  If you give u and v
 970//   as single numbers you'll get a single point back.  If you give u and v as lists or ranges you'll
 971//   get a 2d rectangular array of points.  If one but not both of u and v is a list or range then you'll
 972//   get a list of points.  
 973// Arguments:
 974//   patch = The 2D array of control points for a Bezier patch.
 975//   u = The bezier u parameter (inner list of patch).  Generally between 0 and 1. Can be a list, range or value.
 976//   v = The bezier v parameter (outer list of patch).  Generally between 0 and 1. Can be a list, range or value.
 977// Example(3D):
 978//   patch = [
 979//       [[-50, 50,  0], [-16, 50,  20], [ 16, 50,  20], [50, 50,  0]],
 980//       [[-50, 16, 20], [-16, 16,  40], [ 16, 16,  40], [50, 16, 20]],
 981//       [[-50,-16, 20], [-16,-16,  40], [ 16,-16,  40], [50,-16, 20]],
 982//       [[-50,-50,  0], [-16,-50,  20], [ 16,-50,  20], [50,-50,  0]]
 983//   ];
 984//   debug_bezier_patches(patches=[patch], size=1, showcps=true);
 985//   pt = bezier_patch_points(patch, 0.6, 0.75);
 986//   translate(pt) color("magenta") sphere(d=3, $fn=12);
 987// Example(3D): Getting Multiple Points at Once
 988//   patch = [
 989//       [[-50, 50,  0], [-16, 50,  20], [ 16, 50,  20], [50, 50,  0]],
 990//       [[-50, 16, 20], [-16, 16,  40], [ 16, 16,  40], [50, 16, 20]],
 991//       [[-50,-16, 20], [-16,-16,  40], [ 16,-16,  40], [50,-16, 20]],
 992//       [[-50,-50,  0], [-16,-50,  20], [ 16,-50,  20], [50,-50,  0]]
 993//   ];
 994//   debug_bezier_patches(patches=[patch], size=1, showcps=true);
 995//   pts = bezier_patch_points(patch, [0:0.2:1], [0:0.2:1]);
 996//   for (row=pts) move_copies(row) color("magenta") sphere(d=3, $fn=12);
 997function bezier_patch_points(patch, u, v) =
 998    assert(is_range(u) || is_vector(u) || is_finite(u), "Input u is invalid")
 999    assert(is_range(v) || is_vector(v) || is_finite(v), "Input v is invalid")
1000      !is_num(u) && !is_num(v) ?
1001            let(
1002                vbezes = [for (i = idx(patch[0])) bezier_points(column(patch,i), u)]
1003            )
1004            [for (i = idx(vbezes[0])) bezier_points(column(vbezes,i), v)]
1005    : is_num(u) && is_num(v)? bezier_points([for (bez = patch) bezier_points(bez, v)], u)
1006    : is_num(u) ? bezier_patch_points(patch,force_list(u),v)[0]
1007    :             column(bezier_patch_points(patch,u,force_list(v)),0);
1008
1009
1010  
1011
1012function _bezier_rectangle(patch, splinesteps=16, style="default") =
1013    let(
1014        uvals = lerpn(0,1,splinesteps.x+1),
1015        vvals = lerpn(1,0,splinesteps.y+1),
1016        pts = bezier_patch_points(patch, uvals, vvals)
1017    )
1018    vnf_vertex_array(pts, style=style, reverse=false);
1019
1020
1021// Function: bezier_vnf()
1022// Synopsis: Generates a (probably non-manifold) VNF for one or more bezier surface patches.
1023// SynTags: VNF
1024// Topics: Bezier Patches
1025// See Also: bezier_patch_points(), bezier_patch_flat()
1026// Usage:
1027//   vnf = bezier_vnf(patches, [splinesteps], [style]);
1028// Description:
1029//   Convert a patch or list of patches into the corresponding Bezier surface, representing the
1030//   result as a [VNF structure](vnf.scad).  The `splinesteps` argument specifies the sampling grid of
1031//   the surface for each patch by specifying the number of segments on the borders of the surface.
1032//   It can be a scalar, which gives a uniform grid, or
1033//   it can be [USTEPS, VSTEPS], which gives difference spacing in the U and V parameters. 
1034//   Note that the surface you produce may be disconnected and is not necessarily a valid manifold in OpenSCAD.
1035//   You must also ensure that the patches mate exactly along their edges, or the VNF will be invalid.  
1036// Arguments:
1037//   patches = The bezier patch or list of bezier patches to convert into a vnf.
1038//   splinesteps = Number of segments on the border of the bezier surface.  You can specify [USTEPS,VSTEPS].  Default: 16
1039//   style = The style of subdividing the quads into faces.  Valid options are "default", "alt", "min_edge", "quincunx", "convex" and "concave".  See {{vnf_vertex_array()}}.  Default: "default"
1040// Example(3D):
1041//   patch = [
1042//       // u=0,v=0                                         u=1,v=0
1043//       [[-50,-50,  0], [-16,-50,  20], [ 16,-50, -20], [50,-50,  0]],
1044//       [[-50,-16, 20], [-16,-16,  20], [ 16,-16, -20], [50,-16, 20]],
1045//       [[-50, 16, 20], [-16, 16, -20], [ 16, 16,  20], [50, 16, 20]],
1046//       [[-50, 50,  0], [-16, 50, -20], [ 16, 50,  20], [50, 50,  0]],
1047//       // u=0,v=1                                         u=1,v=1
1048//   ];
1049//   vnf = bezier_vnf(patch, splinesteps=16);
1050//   vnf_polyhedron(vnf);
1051// Example(3D,FlatSpin,VPD=444): Combining multiple patches
1052//   patch = 100*[
1053//       // u=0,v=0                                u=1,v=0
1054//       [[0,  0,0], [1/3,  0,  0], [2/3,  0,  0], [1,  0,0]],
1055//       [[0,1/3,0], [1/3,1/3,1/3], [2/3,1/3,1/3], [1,1/3,0]],
1056//       [[0,2/3,0], [1/3,2/3,1/3], [2/3,2/3,1/3], [1,2/3,0]],
1057//       [[0,  1,0], [1/3,  1,  0], [2/3,  1,  0], [1,  1,0]],
1058//       // u=0,v=1                                u=1,v=1
1059//   ];
1060//   fpatch = bezier_patch_flat([100,100]);
1061//   tpatch = translate([-50,-50,50], patch);
1062//   flatpatch = translate([0,0,50], fpatch);
1063//   vnf = bezier_vnf([
1064//                     tpatch,
1065//                     xrot(90, tpatch),
1066//                     xrot(-90, tpatch),
1067//                     xrot(180, tpatch),
1068//                     yrot(90, flatpatch),
1069//                     yrot(-90, tpatch)]);
1070//   vnf_polyhedron(vnf);
1071// Example(3D):
1072//   patch1 = [
1073//       [[18,18,0], [33,  0,  0], [ 67,  0,  0], [ 82, 18,0]],
1074//       [[ 0,40,0], [ 0,  0,100], [100,  0, 20], [100, 40,0]],
1075//       [[ 0,60,0], [ 0,100,100], [100,100, 20], [100, 60,0]],
1076//       [[18,82,0], [33,100,  0], [ 67,100,  0], [ 82, 82,0]],
1077//   ];
1078//   patch2 = [
1079//       [[18,82,0], [33,100,  0], [ 67,100,  0], [ 82, 82,0]],
1080//       [[ 0,60,0], [ 0,100,-50], [100,100,-50], [100, 60,0]],
1081//       [[ 0,40,0], [ 0,  0,-50], [100,  0,-50], [100, 40,0]],
1082//       [[18,18,0], [33,  0,  0], [ 67,  0,  0], [ 82, 18,0]],
1083//   ];
1084//   vnf = bezier_vnf(patches=[patch1, patch2], splinesteps=16);
1085//   vnf_polyhedron(vnf);
1086// Example(3D): Connecting Patches with asymmetric splinesteps.  Note it is fastest to join all the VNFs at once, which happens in vnf_polyhedron, rather than generating intermediate joined partial surfaces.  
1087//   steps = 8;
1088//   edge_patch = [
1089//       // u=0, v=0                    u=1,v=0
1090//       [[-60, 0,-40], [0, 0,-40], [60, 0,-40]],
1091//       [[-60, 0,  0], [0, 0,  0], [60, 0,  0]],
1092//       [[-60,40,  0], [0,40,  0], [60,40,  0]],
1093//       // u=0, v=1                    u=1,v=1
1094//   ];
1095//   corner_patch = [
1096//       // u=0, v=0                    u=1,v=0
1097//       [[ 0, 40,-40], [ 0,  0,-40], [40,  0,-40]],
1098//       [[ 0, 40,  0], [ 0,  0,  0], [40,  0,  0]],
1099//       [[40, 40,  0], [40, 40,  0], [40, 40,  0]],
1100//       // u=0, v=1                    u=1,v=1
1101//   ];
1102//   face_patch = bezier_patch_flat([120,120],orient=LEFT);
1103//   edges = [
1104//       for (axrot=[[0,0,0],[0,90,0],[0,0,90]], xang=[-90:90:180])
1105//           bezier_vnf(
1106//               splinesteps=[steps,1],
1107//               rot(a=axrot,
1108//                   p=rot(a=[xang,0,0],
1109//                       p=translate(v=[0,-100,100],p=edge_patch)
1110//                   )
1111//               )
1112//           )
1113//   ];
1114//   corners = [
1115//       for (xang=[0,180], zang=[-90:90:180])
1116//           bezier_vnf(
1117//               splinesteps=steps,
1118//               rot(a=[xang,0,zang],
1119//                   p=translate(v=[-100,-100,100],p=corner_patch)
1120//               )
1121//           )
1122//   ];
1123//   faces = [
1124//       for (axrot=[[0,0,0],[0,90,0],[0,0,90]], zang=[0,180])
1125//           bezier_vnf(
1126//               splinesteps=1,
1127//               rot(a=axrot,
1128//                   p=zrot(zang,move([-100,0,0], face_patch))
1129//               )
1130//           )
1131//   ];
1132//   vnf_polyhedron(concat(edges,corners,faces));
1133function bezier_vnf(patches=[], splinesteps=16, style="default") =
1134    assert(is_num(splinesteps) || is_vector(splinesteps,2))
1135    assert(all_positive(splinesteps))
1136    let(splinesteps = force_list(splinesteps,2))
1137    is_bezier_patch(patches)? _bezier_rectangle(patches, splinesteps=splinesteps,style=style)
1138  : assert(is_list(patches),"Invalid patch list")
1139    vnf_join(
1140      [
1141        for (patch=patches)
1142          is_bezier_patch(patch)? _bezier_rectangle(patch, splinesteps=splinesteps,style=style)
1143        : assert(false,"Invalid patch list")
1144      ]
1145    );
1146          
1147
1148
1149// Function: bezier_vnf_degenerate_patch()
1150// Synopsis: Generates a VNF for a degenerate bezier surface patch.
1151// SynTags: VNF
1152// Topics: Bezier Patches
1153// See Also: bezier_patch_points(), bezier_patch_flat(), bezier_vnf()
1154// Usage:
1155//   vnf = bezier_vnf_degenerate_patch(patch, [splinesteps], [reverse]);
1156//   vnf_edges = bezier_vnf_degenerate_patch(patch, [splinesteps], [reverse], return_edges=true);
1157// Description:
1158//   Returns a VNF for a degenerate rectangular bezier patch where some of the corners of the patch are
1159//   equal.  If the resulting patch has no faces then returns an empty VNF.  Note that due to the degeneracy,
1160//   the shape of the surface can be triangular even though the underlying patch is a rectangle.  
1161//   If you specify return_edges then the return is a list whose first element is the vnf and whose second
1162//   element lists the edges in the order [left, right, top, bottom], where each list is a list of the actual
1163//   point values, but possibly only a single point if that edge is degenerate.
1164//   The method checks for various types of degeneracy and uses a triangular or partly triangular array of sample points. 
1165//   See examples below for the types of degeneracy detected and how the patch is sampled for those cases.
1166//   Note that splinesteps is the same for both directions of the patch, so it cannot be an array. 
1167// Arguments:
1168//   patch = Patch to process
1169//   splinesteps = Number of segments to produce on each side.  Default: 16
1170//   reverse = reverse direction of faces.  Default: false
1171//   return_edges = if true return the points on the four edges: [left, right, top, bottom].  Default: false
1172// Example(3D,NoAxes): This quartic patch is degenerate at one corner, where a row of control points are equal.  Processing this degenerate patch normally produces excess triangles near the degenerate point. 
1173//   splinesteps=8;
1174//   patch=[
1175//         repeat([-12.5, 12.5, 15],5),
1176//          [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]],
1177//          [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]],
1178//          [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]],
1179//          [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]]
1180//         ];
1181//   vnf_wireframe((bezier_vnf(patch, splinesteps)),width=0.1);
1182//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1183// Example(3D,NoAxes): With bezier_vnf_degenerate_patch the degenerate point does not have excess triangles.  The top half of the patch decreases the number of sampled points by 2 for each row.  
1184//   splinesteps=8;
1185//   patch=[
1186//          repeat([-12.5, 12.5, 15],5),
1187//          [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]],
1188//          [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]],
1189//          [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]],
1190//          [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]]
1191//         ];
1192//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1193//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1194// Example(3D,NoAxes): With splinesteps odd you get one "odd" row where the point count decreases by 1 instead of 2.  You may prefer even values for splinesteps to avoid this. 
1195//   splinesteps=7;
1196//   patch=[
1197//          repeat([-12.5, 12.5, 15],5),
1198//          [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]],
1199//          [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]],
1200//          [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]],
1201//          [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]]
1202//         ];
1203//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1204//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1205// Example(3D,NoAxes): A more extreme degeneracy occurs when the top half of a patch is degenerate to a line.  (For odd length patches the middle row must be degenerate to trigger this style.)  In this case the number of points in each row decreases by 1 for every row.  It doesn't matter if splinesteps is odd or even. 
1206//   splinesteps=8;
1207//   patch = [[[10, 0, 0], [10, -10.4, 0], [10, -20.8, 0], [1.876, -14.30, 0], [-6.24, -7.8, 0]],
1208//            [[5, 0, 0], [5, -5.2, 0], [5, -10.4, 0], [0.938, -7.15, 0], [-3.12, -3.9, 0]],
1209//            repeat([0,0,0],5),
1210//            repeat([0,0,5],5),
1211//            repeat([0,0,10],5)
1212//           ];
1213//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1214//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1215// Example(3D,NoScales): Here is a degenerate cubic patch.
1216//   splinesteps=8;
1217//   patch = [ [ [-20,0,0],  [-10,0,0],[0,10,0],[0,20,0] ],
1218//             [ [-20,0,10], [-10,0,10],[0,10,10],[0,20,10]],
1219//             [ [-10,0,20], [-5,0,20], [0,5,20], [0,10,20]],
1220//              repeat([0,0,30],4)
1221//               ];
1222//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1223//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1224// Example(3D,NoScales): A more extreme degenerate cubic patch, where two rows are equal.
1225//   splinesteps=8;
1226//   patch = [ [ [-20,0,0], [-10,0,0],[0,10,0],[0,20,0] ],
1227//             [ [-20,0,10], [-10,0,10],[0,10,10],[0,20,10] ],
1228//              repeat([-10,10,20],4),
1229//              repeat([-10,10,30],4)          
1230//           ];
1231//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1232//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1233// Example(3D,NoScales): Quadratic patch degenerate at the right side:
1234//   splinesteps=8;
1235//   patch = [[[0, -10, 0],[10, -5, 0],[20, 0, 0]],
1236//            [[0, 0, 0],  [10, 0, 0], [20, 0, 0]],
1237//            [[0, 0, 10], [10, 0, 5], [20, 0, 0]]];
1238//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1239//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1240// Example(3D,NoAxes): Cubic patch degenerate at both ends.  In this case the point count changes by 2 at every row.  
1241//   splinesteps=8;
1242//   patch = [
1243//            repeat([10,-10,0],4),
1244//            [ [-20,0,0], [-1,0,0],[0,10,0],[0,20,0] ],
1245//            [ [-20,0,10], [-10,0,10],[0,10,10],[0,20,10] ],
1246//            repeat([-10,10,20],4),
1247//           ];
1248//   vnf_wireframe(bezier_vnf_degenerate_patch(patch, splinesteps),width=0.1);
1249//   color("red")move_copies(flatten(patch)) sphere(r=0.3,$fn=9);
1250function bezier_vnf_degenerate_patch(patch, splinesteps=16, reverse=false, return_edges=false) =
1251    !return_edges ? bezier_vnf_degenerate_patch(patch, splinesteps, reverse, true)[0] :
1252    assert(is_bezier_patch(patch), "Input is not a Bezier patch")
1253    assert(is_int(splinesteps) && splinesteps>0, "splinesteps must be a positive integer")
1254    let(
1255        row_degen = [for(row=patch) all_equal(row,eps=EPSILON)],
1256        col_degen = [for(col=transpose(patch)) all_equal(col,eps=EPSILON)],
1257        top_degen = row_degen[0],
1258        bot_degen = last(row_degen),
1259        left_degen = col_degen[0],
1260        right_degen = last(col_degen),
1261        samplepts = lerpn(0,1,splinesteps+1)
1262    )
1263    all(row_degen) && all(col_degen) ?  // fully degenerate case
1264        [EMPTY_VNF, repeat([patch[0][0]],4)] :
1265    all(row_degen) ?                         // degenerate to a line (top to bottom)
1266        let(pts = bezier_points(column(patch,0), samplepts))
1267        [EMPTY_VNF, [pts,pts,[pts[0]],[last(pts)]]] :
1268    all(col_degen) ?                         // degenerate to a line (left to right)
1269        let(pts = bezier_points(patch[0], samplepts))
1270        [EMPTY_VNF, [[pts[0]], [last(pts)], pts, pts]] :
1271    !top_degen && !bot_degen && !left_degen && !right_degen ?       // non-degenerate case
1272       let(pts = bezier_patch_points(patch, samplepts, samplepts))
1273       [
1274        vnf_vertex_array(pts, reverse=!reverse),
1275        [column(pts,0), column(pts,len(pts)-1), pts[0], last(pts)]
1276       ] :
1277    top_degen && bot_degen ?
1278       let(
1279            rowcount = [
1280                        each list([3:2:splinesteps]),
1281                        if (splinesteps%2==0) splinesteps+1,
1282                        each reverse(list([3:2:splinesteps]))
1283                       ],
1284            bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(column(patch,i), samplepts)],
1285            pts = [
1286                  [bpatch[0][0]],
1287                  for(j=[0:splinesteps-2]) bezier_points(column(bpatch,j+1), lerpn(0,1,rowcount[j])),
1288                  [last(bpatch[0])]
1289                  ],
1290            vnf = vnf_tri_array(pts, reverse=!reverse)
1291         ) [
1292            vnf,
1293            [
1294             column(pts,0),
1295             [for(row=pts) last(row)],
1296             pts[0],
1297             last(pts),
1298            ]
1299          ]  :    
1300    bot_degen ?                                           // only bottom is degenerate
1301       let(
1302           result = bezier_vnf_degenerate_patch(reverse(patch), splinesteps=splinesteps, reverse=!reverse, return_edges=true)
1303       )
1304       [
1305          result[0],
1306          [reverse(result[1][0]), reverse(result[1][1]), (result[1][3]), (result[1][2])]
1307       ] :
1308    top_degen ?                                          // only top is degenerate
1309       let(
1310           full_degen = len(patch)>=4 && all(select(row_degen,1,ceil(len(patch)/2-1))),
1311           rowmax = full_degen ? count(splinesteps+1) :
1312                                 [for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps],
1313           bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(column(patch,i), samplepts)],
1314           pts = [
1315                  [bpatch[0][0]],
1316                  for(j=[1:splinesteps]) bezier_points(column(bpatch,j), lerpn(0,1,rowmax[j]+1))
1317                 ],
1318           vnf = vnf_tri_array(pts, reverse=!reverse)
1319        ) [
1320            vnf,
1321            [
1322             column(pts,0),
1323             [for(row=pts) last(row)],
1324             pts[0],
1325             last(pts),
1326            ]
1327          ] :
1328      // must have left or right degeneracy, so transpose and recurse
1329      let(
1330          result = bezier_vnf_degenerate_patch(transpose(patch), splinesteps=splinesteps, reverse=!reverse, return_edges=true)
1331      )
1332      [result[0],
1333       select(result[1],[2,3,0,1])
1334      ];
1335
1336
1337// Function: bezier_patch_normals()
1338// Synopsis: Computes surface normals for one or more places on a bezier surface patch.
1339// Topics: Bezier Patches
1340// See Also: bezier_patch_points(), bezier_points(), bezier_curve(), bezpath_curve()
1341// Usage:
1342//   n = bezier_patch_normals(patch, u, v);
1343//   ngrid = bezier_patch_normals(patch, LIST, LIST);
1344//   ngrid = bezier_patch_normals(patch, RANGE, RANGE);
1345// Description:
1346//   Compute the unit normal vector to a bezier patch at the listed point set.  The bezier patch must be a rectangular array of
1347//   points, and the normal will be computed at all the (u,v) pairs that you specify.  If you give u and v
1348//   as single numbers you'll get a single point back.  If you give u and v as lists or ranges you'll
1349//   get a 2d rectangular array of points.  If one but not both of u and v is a list or range then you'll
1350//   get a list of points.
1351//   .
1352//   This function works by computing the cross product of the tangents.  In some degenerate cases the one of the tangents
1353//   can be zero, so the normal vector does not exist.  In this case, undef is returned.  Another degenerate case
1354//   occurs when the tangents are parallel, or nearly parallel.  In this case you will get a unit vector returned but it will not
1355//   be the correct normal vector.   This can happen if you use a degenerate patch, or if you give two of the edges of your patch a smooth "corner"
1356//   so that the u and v directions are parallel at the corner.  
1357// Arguments:
1358//   patch = The 2D array of control points for a Bezier patch.
1359//   u = The bezier u parameter (inner list of patch).  Generally between 0 and 1. Can be a list, range or value.
1360//   v = The bezier v parameter (outer list of patch).  Generally between 0 and 1. Can be a list, range or value.
1361// Example(3D,Med,VPR=[71.1,0,155.9],VPD=292.705,VPT=[20.4724,38.7273,22.7683],NoAxes): Normal vectors on a patch
1362//   patch = [
1363//        // u=0,v=0                                         u=1,v=0
1364//        [[-50,-50,  0], [-16,-50,  20], [ 16,-50, -20], [50,-50,  0]],
1365//        [[-50,-16, 40], [-16,-16,  20], [ 16,-16, -20], [50,-16, 70]],
1366//        [[-50, 16, 20], [-16, 16, -20], [ 16, 37,  20], [70, 16, 20]],
1367//        [[-50, 50,  0], [73, 50, -40], [ 16, 50,  20], [50, 50,  0]],
1368//        // u=0,v=1                                         u=1,v=1
1369//   ];
1370//   vnf_polyhedron(bezier_vnf(patch,splinesteps=30));
1371//   uv = lerpn(0,1,12);
1372//   pts = bezier_patch_points(patch, uv, uv);
1373//   normals = bezier_patch_normals(patch, uv, uv);
1374//     for(i=idx(uv),j=idx(uv)){
1375//        stroke([pts[i][j],pts[i][j]-6*normals[i][j]], width=0.5,
1376//               endcap1="dot",endcap2="arrow2",color="blue");
1377//   }
1378// Example(3D,NoAxes,Med,VPR=[72.5,0,288.9],VPD=192.044,VPT=[51.6089,48.118,5.89088]): This example gives invalid normal vectors at the four corners of the patch where the u and v directions are parallel.  You can see how the line of triangulation is approaching parallel at the edge, and the invalid vectors in red point in a completely incorrect direction.  
1379//   patch = [
1380//       [[18,18,0], [33,  0,  0], [ 67,  0,  0], [ 82, 18,0]],
1381//       [[ 0,40,0], [ 0,  0,100], [100,  0, 20], [100, 40,0]],
1382//       [[ 0,60,0], [ 0,100,100], [100,100, 20], [100, 60,0]],
1383//       [[18,82,0], [33,100,  0], [ 67,100,  0], [ 82, 82,0]],
1384//   ];
1385//   vnf_polyhedron(bezier_vnf(patch,splinesteps=30));
1386//   uv = lerpn(0,1,7);
1387//   pts = bezier_patch_points(patch, uv, uv);
1388//   normals = bezier_patch_normals(patch, uv, uv);
1389//     for(i=idx(uv),j=idx(uv)){
1390//        color=((uv[i]==0 || uv[i]==1) && (uv[j]==0 || uv[j]==1))
1391//             ? "red" : "blue";
1392//        stroke([pts[i][j],pts[i][j]-8*normals[i][j]], width=0.5,
1393//               endcap1="dot",endcap2="arrow2",color=color);
1394//   }
1395// Example(3D,Med,NoAxes,VPR=[56.4,0,71.9],VPD=66.9616,VPT=[10.2954,1.33721,19.4484]): This degenerate patch has normals everywhere, but computation of the normal fails at the point of degeneracy, the top corner.  
1396//    patch=[
1397//             repeat([-12.5, 12.5, 15],5),
1398//              [[-6.25, 11.25, 15], [-6.25, 8.75, 15], [-6.25, 6.25, 15], [-8.75, 6.25, 15], [-11.25, 6.25, 15]],
1399//              [[0, 10, 15], [0, 5, 15], [0, 0, 15], [-5, 0, 15], [-10, 0, 15]],
1400//              [[0, 10, 8.75], [0, 5, 8.75], [0, 0, 8.75], [-5, 0, 8.75], [-10, 0, 8.75]],
1401//              [[0, 10, 2.5], [0, 5, 2.5], [0, 0, 2.5], [-5, 0, 2.5], [-10, 0, 2.5]]
1402//             ];
1403//    vnf_polyhedron(bezier_vnf(patch, 32));
1404//    uv = lerpn(0,1,8);
1405//    pts = bezier_patch_points(patch, uv, uv);
1406//    normals = bezier_patch_normals(patch, uv, uv);
1407//      for(i=idx(uv),j=idx(uv)){
1408//        if (is_def(normals[i][j]))
1409//          stroke([pts[i][j],pts[i][j]-2*normals[i][j]], width=0.1,
1410//                 endcap1="dot",endcap2="arrow2",color="blue");
1411//    }
1412// Example(3D,Med,NoAxes,VPR=[48,0,23.6],VPD=32.0275,VPT=[-0.145727,-0.0532125,1.74224]): This example has a singularities where the tangent lines don't exist, so the normal will be undef at those points.  
1413//    pts1 = [ [-5,0,0], [5,0,5], [-5,0,5], [5,0,0] ];
1414//    pts2 = [ [0,-5,0], [0,5,5], [0,-5,5], [0,5,0] ];
1415//    patch = [for(i=[0:3])
1416//            [for(j=[0:3]) pts1[i]+pts2[j] ] ];
1417//    vnf_polyhedron(bezier_vnf(patch, 163));
1418//    uv = [0,.1,.2,.3,.7,.8,.9,1];//lerpn(0,1,8);
1419//    pts = bezier_patch_points(patch, uv, uv);
1420//    normals = bezier_patch_normals(patch, uv, uv);
1421//    for(i=idx(uv),j=idx(uv))
1422//      stroke([pts[i][j],pts[i][j]+2*normals[i][j]], width=0.08,
1423//               endcap1="dot",endcap2="arrow2",color="blue");
1424  
1425function bezier_patch_normals(patch, u, v) =
1426    assert(is_range(u) || is_vector(u) || is_finite(u), "Input u is invalid")
1427    assert(is_range(v) || is_vector(v) || is_finite(v), "Input v is invalid")
1428      !is_num(u) && !is_num(v) ?
1429          let(
1430              vbezes = [for (i = idx(patch[0])) bezier_points(column(patch,i), u)],
1431              dvbezes = [for (i = idx(patch[0])) bezier_derivative(column(patch,i), u)],
1432              v_tangent = [for (i = idx(vbezes[0])) bezier_derivative(column(vbezes,i), v)],
1433              u_tangent = [for (i = idx(vbezes[0])) bezier_points(column(dvbezes,i), v)]
1434          )
1435          [for(i=idx(u_tangent)) [for(j=idx(u_tangent[0])) unit(cross(u_tangent[i][j],v_tangent[i][j]),undef)]]
1436    : is_num(u) && is_num(v)?
1437          let(
1438                du = bezier_derivative([for (bez = patch) bezier_points(bez, v)], u),
1439                dv = bezier_points([for (bez = patch) bezier_derivative(bez, v)], u)
1440          )
1441          unit(cross(du,dv),undef)
1442    : is_num(u) ? bezier_patch_normals(patch,force_list(u),v)[0]
1443    :             column(bezier_patch_normals(patch,u,force_list(v)),0);
1444
1445
1446
1447// Section: Debugging Beziers
1448
1449
1450// Module: debug_bezier()
1451// Synopsis: Shows a bezier path and its associated control points.
1452// SynTags: Geom
1453// Topics: Bezier Paths, Debugging
1454// See Also: bezpath_curve()
1455// Usage:
1456//   debug_bezier(bez, [size], [N=]);
1457// Description:
1458//   Renders 2D or 3D bezier paths and their associated control points.
1459//   Useful for debugging bezier paths.
1460// Arguments:
1461//   bez = the array of points in the bezier.
1462//   size = diameter of the lines drawn.
1463//   ---
1464//   N = Mark the first and every Nth vertex after in a different color and shape.
1465// Example(2D):
1466//   bez = [
1467//       [-10,   0],  [-15,  -5],
1468//       [ -5, -10],  [  0, -10],  [ 5, -10],
1469//       [ 14,  -5],  [ 15,   0],  [16,   5],
1470//       [  5,  10],  [  0,  10]
1471//   ];
1472//   debug_bezier(bez, N=3, width=0.5);
1473module debug_bezier(bezpath, width=1, N=3) {
1474    no_children($children);
1475    check = 
1476      assert(is_path(bezpath),"bezpath must be a path")
1477      assert(is_int(N) && N>0, "N must be a positive integer")
1478      assert(len(bezpath)%N == 1, str("A degree ",N," bezier path shound have a multiple of ",N," points in it, plus 1."));
1479    $fn=8;
1480    stroke(bezpath_curve(bezpath, N=N), width=width, color="cyan");
1481    color("green")
1482      if (N!=3) 
1483           stroke(bezpath, width=width);
1484      else 
1485           for(i=[1:3:len(bezpath)]) stroke(select(bezpath,max(0,i-2), min(len(bezpath)-1,i)), width=width);
1486    twodim = len(bezpath[0])==2;
1487    color("red") move_copies(bezpath)
1488      if ($idx % N !=0)
1489          if (twodim){
1490            rect([width/2, width*3]);
1491            rect([width*3, width/2]);
1492          } else {
1493           zcyl(d=width/2, h=width*3);
1494           xcyl(d=width/2, h=width*3);
1495           ycyl(d=width/2, h=width*3);
1496        }
1497    color("blue") move_copies(bezpath)
1498      if ($idx % N ==0)
1499        if (twodim) circle(d=width*2.25); else sphere(d=width*2.25);
1500    if (twodim) color("red") move_copies(bezpath)
1501      if ($idx % N !=0) circle(d=width/2);
1502}
1503
1504
1505// Module: debug_bezier_patches()
1506// Synopsis: Shows a bezier surface patch and its associated control points.
1507// SynTags: Geom
1508// Topics: Bezier Patches, Debugging
1509// See Also: bezier_patch_points(), bezier_patch_flat(), bezier_vnf()
1510// Usage:
1511//   debug_bezier_patches(patches, [size=], [splinesteps=], [showcps=], [showdots=], [showpatch=], [convexity=], [style=]);
1512// Description:
1513//   Shows the surface, and optionally, control points of a list of bezier patches.
1514// Arguments:
1515//   patches = A list of rectangular bezier patches.
1516//   ---
1517//   splinesteps = Number of segments to divide each bezier curve into. Default: 16
1518//   showcps = If true, show the controlpoints as well as the surface.  Default: true.
1519//   showdots = If true, shows the calculated surface vertices.  Default: false.
1520//   showpatch = If true, shows the surface faces.  Default: true.
1521//   size = Size to show control points and lines.  Default: 1% of the maximum side length of a box bounding the patch.
1522//   style = The style of subdividing the quads into faces.  Valid options are "default", "alt", and "quincunx".
1523//   convexity = Max number of times a line could intersect a wall of the shape.
1524// Example:
1525//   patch1 = [
1526//       [[15,15,0], [33,  0,  0], [ 67,  0,  0], [ 85, 15,0]],
1527//       [[ 0,33,0], [33, 33, 50], [ 67, 33, 50], [100, 33,0]],
1528//       [[ 0,67,0], [33, 67, 50], [ 67, 67, 50], [100, 67,0]],
1529//       [[15,85,0], [33,100,  0], [ 67,100,  0], [ 85, 85,0]],
1530//   ];
1531//   patch2 = [
1532//       [[15,85,0], [33,100,  0], [ 67,100,  0], [ 85, 85,0]],
1533//       [[ 0,67,0], [33, 67,-50], [ 67, 67,-50], [100, 67,0]],
1534//       [[ 0,33,0], [33, 33,-50], [ 67, 33,-50], [100, 33,0]],
1535//       [[15,15,0], [33,  0,  0], [ 67,  0,  0], [ 85, 15,0]],
1536//   ];
1537//   debug_bezier_patches(patches=[patch1, patch2], splinesteps=8, showcps=true);
1538module debug_bezier_patches(patches=[], size, splinesteps=16, showcps=true, showdots=false, showpatch=true, convexity=10, style="default")
1539{
1540    no_children($children);
1541    assert(is_undef(size)||is_num(size));
1542    assert(is_int(splinesteps) && splinesteps>0);
1543    assert(is_list(patches) && all([for (patch=patches) is_bezier_patch(patch)]));
1544    assert(is_bool(showcps));
1545    assert(is_bool(showdots));
1546    assert(is_bool(showpatch));
1547    assert(is_int(convexity) && convexity>0);
1548    for (patch = patches) {
1549        size = is_num(size)? size :
1550               let( bounds = pointlist_bounds(flatten(patch)) )
1551               max(bounds[1]-bounds[0])*0.01;
1552        if (showcps) {
1553            move_copies(flatten(patch)) color("red") sphere(d=size*2);
1554            color("cyan") 
1555                for (i=[0:1:len(patch)-1], j=[0:1:len(patch[i])-1]) {
1556                        if (i<len(patch)-1) extrude_from_to(patch[i][j], patch[i+1][j]) circle(d=size);
1557                        if (j<len(patch[i])-1) extrude_from_to(patch[i][j], patch[i][j+1]) circle(d=size);
1558                }        
1559        }
1560        if (showpatch || showdots){
1561            vnf = bezier_vnf(patch, splinesteps=splinesteps, style=style);
1562            if (showpatch) vnf_polyhedron(vnf, convexity=convexity);
1563            if (showdots) color("blue") move_copies(vnf[0]) sphere(d=size);
1564        }
1565    }
1566}
1567
1568
1569// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap