1 /* bezier.i 2 * compute points on Bezier curve 3 */ 4 /* Copyright (c) 2013, David H. Munro. 5 * All rights reserved. 6 * This file is part of yorick (http://yorick.github.com). 7 * Read the accompanying LICENSE file for details. 8 */ 9 10 func bezier(p, s, n, w=) 11 /* DOCUMENT bezier(p, s) 12 * or bezier(p, s, n) 13 * 14 * Return points on a Bezier curve defined by the array of control 15 * points P, at parameter values S. The Bezier curve is of order N, 16 * which defaults to 3. P is a d-by-(m*N+1) dimensional array 17 * of control points in d-dimensions. (That is, d=2 is a planar curve, 18 * d=3 is a curve in 3D space, and so on.) The points 1, 1+N, 1+2*N, 19 * 1+3*N, etc. are on the curve, while the N-1 intermediate points 20 * are Bezier control points. The array S may have any dimensionality; 21 * the return value will have dimensions d-by-dimsof(S). The values 22 * of S must lie between 0 and m, where S=0 represents P(,1), S=1 23 * represents P(,2), S=2 represents P(,3), and so on, and intermediate 24 * values of S are on the Bezier curve connecting those integer points 25 * defined by the intermediate control points. 26 * 27 * As a shorthand, if S is an integer scalar greater than 1, bezier 28 * gives you S equally spaced points from 0 to m, S = span(0,m,S). 29 * 30 * If the endpoints of one interval are P[0] and P[N], and the control 31 * points are P[1], P[2], ..., P[N-1], then the Bezier curve of order N 32 * is defined by: 33 * P(s) = sum[k=0,N]{ P[k] * (N k)*(1-s)^(N-k)*s^k } 34 * where (N k) = N!/((N-1)!*k!) is the binomial coefficient, where s 35 * goes from 0 to 1 takes P(s) from P[0] to P[N] along the Bezier curve 36 * of order N. 37 * 38 * The default case n=3 for d=1 produces a piecewise cubic curve, which 39 * you can also produce using the spline or splinef function. 40 * 41 * The w= keyword can be used to supply weights W to produce a rational 42 * Bezier curve with control points P and weights W: 43 * q = bezier(P, s, w=W); 44 * is equivalent to: 45 * pw = W(-:1:dimsof(P)(2)+1,); 46 * pw(1:-1,..) *= P; 47 * q = bezier(pw, s); 48 * q = q(1:-1,..) / q(0:0,..); 49 * Rational Bezier curves of order 2 exactly represent ellipses, 50 * parabolas, or hyperbolas, when the weights w are 1 for points on the 51 * curve and w<1, w=1, or w>1, respectively, for intermediate points. 52 * For example, bezier([[0,1],[1,1],[1,0]],w=[1,sqrt(.5),1], 100, 2) is 53 * 100 points along the first quadrant of the unit circle. 54 * 55 * SEE ALSO: span, spanarc, spline, splinef, interp 56 */ 57 { 58 if (is_void(n)) n = 3; 59 d = dimsof(p); 60 m = (numberof(d) && d(1)==2)? (d(3)-1)/n : -1; 61 if (m<=0 || d(3)!=m*n+1) 62 error, "wrong number of control points, need npts%order = 1"; 63 if (!is_void(w)) { 64 if (numberof(w)!=d(3) || dimsof(w)(1)!=1) 65 error, "wrong number of point weights w= for rational bezier"; 66 w = w(-:1:d(2)+1,..); 67 w(1:-1,) *= p; 68 p = w; 69 w = 1; 70 } 71 if (numberof(s)==1 && structof(s+0)==long) s = span(0, m, s); 72 i = min(long(max(s,0.)), m-1); 73 s = (s-i)(-,-,..); 74 /* iterated linear interpolation, modified De Casteljau's algorithm */ 75 p = p(, indgen(n+1) + n*i(-,..)); 76 for (k=0 ; k<n ; ++k) p = p(,1:-1,..) + p(,dif,..) * s; 77 p = p(,1,..); 78 if (w) p = p(1:-1,..) / p(0:0,..); 79 return p; 80 } 81