FE Verification using photoelasticity and constrained shape optimization

NURBS_toolbox.zip

Home >Engineering homework help >Mechanical Engineering homework help >FE Verification using photoelasticity and constrained shape optimization

nurbs_toolbox/veccross.m

function cross = veccross(vec1,vec2) % % Function Name: % % veccross - The cross product of two vectors. % % Calling Sequence: % % cross = veccross(vec1,vec2); % % Parameters: % % vec1 : An array of column vectors represented by a matrix of % vec2 size (dim,nv), where is the dimension of the vector and % nv the number of vectors. % % cross : Array of column vectors, each element is corresponding % to the cross product of the respective components in vec1 % and vec2. % % Description: % % Cross product of two vectors. % % Examples: % % Determine the cross products of: % (2.3,3.4,5.6) and (1.2,4.5,1.2) % (5.1,0.0,2.3) and (2.5,3.2,4.0) % % cross = veccross([2.3 5.1; 3.4 0.0; 5.6 2.3],[1.2 2.5; 4.5 3.2; 1.2 4.0]); % D.M. Spink % Copyright (c) 2000. if size(vec1,1) == 2 % 2D vector cross = zeros(size(vec1)); cross(3,:) = vec1(1,:).*vec2(2,:)-vec1(2,:).*vec2(1,:); else % 3D vector cross = [vec1(2,:).*vec2(3,:)-vec1(3,:).*vec2(2,:); vec1(3,:).*vec2(1,:)-vec1(1,:).*vec2(3,:); vec1(1,:).*vec2(2,:)-vec1(2,:).*vec2(1,:)]; end

nurbs_toolbox/vecdot.m

function dot = vecdot(vec1,vec2) % % Function Name: % % vecdot - The dot product of two vectors. % % Calling Sequence: % % dot = vecdot(vec1,vec2); % % Parameters: % % vec1 : An array of column vectors represented by a matrix of % vec2 size (dim,nv), where is the dimension of the vector and % nv the number of vectors. % % dot : Row vector of scalars, each element corresponding to % the dot product of the respective components in vec1 and % vec2. % % Description: % % Scalar dot product of two vectors. % % Examples: % % Determine the dot product of % (2.3,3.4,5.6) and (1.2,4.5,1.2) % (5.1,0.0,2.3) and (2.5,3.2,4.0) % % dot = vecdot([2.3 5.1; 3.4 0.0; 5.6 2.3],[1.2 2.5; 4.5 3.2; 1.2 4.0]); % D.M. Spink % Copyright (c) 2000 dot = sum(vec1.*vec2);

nurbs_toolbox/vecmag2.m

function mag = vecmag2(vec) % % Function Name: % % vecmag - Squared magnitude of the vectors % % Calling Sequence: % % mvec = vecmag2(vec) % % Parameters: % % vec : An array of column vectors represented by a matrix of % size (dim,nv), where is the dimension of the vector and % nv the number of vectors. % % mvec : Squared magnitude of the vectors, vector of size (1,nv). % % Description: % % Determines the squared magnitude of the vectors. % % Examples: % % Find the squared magnitude of the two vectors (0.0,2.0,1.3) % and (1.5,3.4,2.3) % % mvec = vecmag2([0.0 1.5; 2.0 3.4; 1.3 2.3]); % D.M. Spink % Copyright (c) 2000. mag = sum(vec.^2);

nurbs_toolbox/vecmag.m

function mag = vecmag(vec) % % Function Name: % % vecmag - Magnitude of the vectors % % Calling Sequence: % % mvec = vecmag(vec) % % Parameters: % % vec : An array of column vectors represented by a matrix of % size (dim,nv), where is the dimension of the vector and % nv the number of vectors. % % mvec : Magnitude of the vectors, vector of size (1,nv). % % Description: % % Determines the magnitude of the vectors. % % Examples: % % Find the magnitude of the two vectors (0.0,2.0,1.3) and (1.5,3.4,2.3) % % mvec = vecmag([0.0 1.5; 2.0 3.4; 1.3 2.3]); % D.M. Spink % Copyright (c) 2000. mag = sqrt(sum(vec.^2));

nurbs_toolbox/vecnorm.m

function nvec = vecnorm(vec) % % Function Name: % % vecnorm - Normalise the vectors. % % Calling Sequence: % % nvec = vecnorn(vec); % % Parameters: % % vec : An array of column vectors represented by a matrix of % size (dim,nv), where is the dimension of the vector and % nv the number of vectors. % % nvec : Normalised vectors, matrix the smae size as vec. % % Description: % % Normalises the array of vectors, returning the unit vectors. % % Examples: % % Normalise the two vectors (0.0,2.0,1.3) and (1.5,3.4,2.3) % % nvec = vecnorm([0.0 1.5; 2.0 3.4; 1.3 2.3]); % D.M. Spink % Copyright (c) 2000. nvec = vec./repmat(sqrt(sum(vec.^2)),[size(vec,1) ones(1,ndims(vec)-1)]);

nurbs_toolbox/vecrotx.m

function rx = vecrotx(angle) % % Function Name: % % vecrotx - Transformation matrix for a rotation around the x axis. % % Calling Sequence: % % rx = vecrotx(angle); % % Parameters: % % angle : rotation angle defined in radians % % rx : (4x4) Transformation matrix. % % % Description: % % Return the (4x4) Transformation matrix for a rotation about the x axis % by the defined angle. % % The matrix is: % % [ 1 0 0 0] % [ 0 cos(angle) -sin(angle) 0] % [ 0 sin(angle) cos(angle) 0] % [ 0 0 0 1] % % Examples: % % Rotate the NURBS line (0.0 0.0 0.0) - (3.0 3.0 3.0) by 45 degrees % around the x-axis % % line = nrbline([0.0 0.0 0.0],[3.0 3.0 3.0]); % trans = vecrotx(%pi/4); % rline = nrbtform(line, trans); % % See: % % nrbtform % Dr D.M. Spink % Copyright (c) 2000. sn = sin(angle); cn = cos(angle); rx = [1 0 0 0; 0 cn -sn 0; 0 sn cn 0; 0 0 0 1];

nurbs_toolbox/vecroty.m

function ry = vecroty(angle) % % Function Name: % % vecroty - Transformation matrix for a rotation around the y axis. % % Calling Sequence: % % ry = vecroty(angle); % % Parameters: % % angle : rotation angle defined in radians % % ry : (4x4) Transformation matrix. % % % Description: % % Return the (4x4) Transformation matrix for a rotation about the y axis % by the defined angle. % % The matrix is: % % [ cos(angle) 0 sin(angle) 0] % [ 0 1 0 0] % [ -sin(angle) 0 cos(angle) 0] % [ 0 0 0 1] % % Examples: % % Rotate the NURBS line (0.0 0.0 0.0) - (3.0 3.0 3.0) by 45 degrees % around the y-axis % % line = nrbline([0.0 0.0 0.0],[3.0 3.0 3.0]); % trans = vecroty(%pi/4); % rline = nrbtform(line, trans); % % See: % % nrbtform % Dr D.M. Spink % Copyright (c) 2000. sn = sin(angle); cn = cos(angle); ry = [cn 0 sn 0; 0 1 0 0; -sn 0 cn 0; 0 0 0 1];

nurbs_toolbox/vecrotz.m

function rz = vecrotz(angle) % % Function Name: % % vecrotz - Transformation matrix for a rotation around the z axis. % % Calling Sequence: % % rz = vecrotz(angle); % % Parameters: % % angle : rotation angle defined in radians % % rz : (4x4) Transformation matrix. % % % Description: % % Return the (4x4) Transformation matrix for a rotation about the z axis % by the defined angle. % % The matrix is: % % [ cos(angle) -sin(angle) 0 0] % [ -sin(angle) cos(angle) 0 0] % [ 0 0 1 0] % [ 0 0 0 1] % % Examples: % % Rotate the NURBS line (0.0 0.0 0.0) - (3.0 3.0 3.0) by 45 degrees % around the z-axis % % line = nrbline([0.0 0.0 0.0],[3.0 3.0 3.0]); % trans = vecrotz(%pi/4); % rline = nrbtform(line, trans); % % See: % % nrbtform % Dr D.M. Spink % Copyright (c) 2000. sn = sin(angle); cn = cos(angle); rz = [cn -sn 0 0; sn cn 0 0; 0 0 1 0; 0 0 0 1];

nurbs_toolbox/vecscale.m

function ss = vecscale(vector) % % Function Name: % % vecscale - Transformation matrix for a scaling. % % Calling Sequence: % % ss = vecscale(svec) % % Parameters: % % svec : A vectors defining the scaling along the x,y and z axes. % i.e. [sx, sy, sy] % % ss : Scaling Transformation Matrix % % Description: % % Returns a (4x4) Transformation matrix for scaling. % % The matrix is: % % [ sx 0 0 0] % [ 0 sy 0 0] % [ 0 0 sz 0] % [ 0 0 0 1] % % Example: % % Scale up the NURBS line (0.0,0.0,0.0) - (1.0,1.0,1.0) by 3 along % the x-axis, 2 along the y-axis and 4 along the z-axis. % % line = nrbline([0.0 0.0 0.0],[1.0 1.0 1.0]); % trans = vecscale([3.0 2.0 4.0]); % sline = nrbtform(line, trans); % % See: % % nrbtform % Dr D.M. Spink % Copyright (c) 2000. if nargin < 1 error('Scaling vector not specified'); end s = [vector(:);0;0]; ss = [s(1) 0 0 0; 0 s(2) 0 0; 0 0 s(3) 0; 0 0 0 1];

nurbs_toolbox/vectrans.m

function dd = vectrans(vector) % % Function Name: % % vectrans - Transformation matrix for a translation. % % Calling Sequence: % % st = vectrans(tvec) % % Parameters: % % tvec : A vectors defining the translation along the x,y and % z axes. i.e. [tx, ty, ty] % % st : Translation Transformation Matrix % % Description: % % Returns a (4x4) Transformation matrix for translation. % % The matrix is: % % [ 1 0 0 tx ] % [ 0 0 0 ty ] % [ 0 0 0 tz ] % [ 0 0 0 1 ] % % Examples: % % Translate the NURBS line (0.0,0.0,0.0) - (1.0,1.0,1.0) by 3 along % the x-axis, 2 along the y-axis and 4 along the z-axis. % % line = nrbline([0.0 0.0 0.0],[1.0 1.0 1.0]); % trans = vectrans([3.0 2.0 4.0]); % tline = nrbtform(line, trans); % % See: % % nrbtform % Dr D.M. Spink % Copyright (c) 2000. if nargin < 1 error('Translation vector required'); end v = [vector(:);0;0]; dd = [1 0 0 v(1); 0 1 0 v(2); 0 0 1 v(3); 0 0 0 1];

nurbs_toolbox/basisfun.m

function N = basisfun(i,u,p,U) % BASISFUN Basis function for B-Spline % ------------------------------------------------------------------------- % ADAPTATION of BASISFUN from C Routine % ------------------------------------------------------------------------- % % Calling Sequence: % % N = basisfun(i,u,p,U) % % INPUT: % % i - knot span ( from FindSpan() ) % u - parametric point % p - spline degree % U - knot sequence % % OUTPUT: % % N - Basis functions vector[p+1] % % Algorithm A2.2 from 'The NURBS BOOK' pg70. % void basisfun(int i, double u, int p, double *U, double *N) { % int j,r; % double saved, temp; i = i + 1; % // work space left = zeros(p+1,1); % double *left = (double*) mxMalloc((p+1)*sizeof(double)); right = zeros(p+1,1); % double *right = (double*) mxMalloc((p+1)*sizeof(double)); N(1) = 1; % N[0] = 1.0; for j=1:p % for (j = 1; j <= p; j++) { left(j+1) = u - U(i+1-j); % left[j] = u - U[i+1-j]; right(j+1) = U(i+j) - u; % right[j] = U[i+j] - u; saved = 0; % saved = 0.0; for r=0:j-1 % for (r = 0; r < j; r++) { temp = N(r+1)/(right(r+2) + left(j-r+1)); % temp = N[r] / (right[r+1] + left[j-r]); N(r+1) = saved + right(r+2)*temp; % N[r] = saved + right[r+1] * temp; saved = left(j-r+1)*temp; % saved = left[j-r] * temp; end % } N(j+1) = saved; % N[j] = saved; end % } % mxFree(left); % mxFree(right); % }

nurbs_toolbox/bspdegelev.m

function [ic,ik] = bspdegelev(d,c,k,t) % % Function Name: % % bspdegevel - Degree elevate a univariate B-Spline. % % Calling Sequence: % % [ic,ik] = bspdegelev(d,c,k,t) % % Parameters: % % d : Degree of the B-Spline. % % c : Control points, matrix of size (dim,nc). % % k : Knot sequence, row vector of size nk. % % t : Raise the B-Spline degree t times. % % ic : Control points of the new B-Spline. % % ik : Knot vector of the new B-Spline. % % Description: % % Degree elevate a univariate B-Spline. This function provides an % interface to a toolbox 'C' routine. [mc,nc] = size(c); % % int bspdegelev(int d, double *c, int mc, int nc, double *k, int nk, % int t, int *nh, double *ic, double *ik) % { % int row,col; % % int ierr = 0; % int i, j, q, s, m, ph, ph2, mpi, mh, r, a, b, cind, oldr, mul; % int n, lbz, rbz, save, tr, kj, first, kind, last, bet, ii; % double inv, ua, ub, numer, den, alf, gam; % double **bezalfs, **bpts, **ebpts, **Nextbpts, *alfs; % % double **ctrl = vec2mat(c, mc, nc); % ic = zeros(mc,nc*(t)); % double **ictrl = vec2mat(ic, mc, nc*(t+1)); % n = nc - 1; % n = nc - 1; % bezalfs = zeros(d+1,d+t+1); % bezalfs = matrix(d+1,d+t+1); bpts = zeros(mc,d+1); % bpts = matrix(mc,d+1); ebpts = zeros(mc,d+t+1); % ebpts = matrix(mc,d+t+1); Nextbpts = zeros(mc,d+1); % Nextbpts = matrix(mc,d+1); alfs = zeros(d,1); % alfs = (double *) mxMalloc(d*sizeof(double)); % m = n + d + 1; % m = n + d + 1; ph = d + t; % ph = d + t; ph2 = floor(ph / 2); % ph2 = ph / 2; % % // compute bezier degree elevation coefficeients bezalfs(1,1) = 1; % bezalfs[0][0] = bezalfs[ph][d] = 1.0; bezalfs(d+1,ph+1) = 1; % for i=1:ph2 % for (i = 1; i <= ph2; i++) { inv = 1/bincoeff(ph,i); % inv = 1.0 / bincoeff(ph,i); mpi = min(d,i); % mpi = min(d,i); % for j=max(0,i-t):mpi % for (j = max(0,i-t); j <= mpi; j++) bezalfs(j+1,i+1) = inv*bincoeff(d,j)*bincoeff(t,i-j); % bezalfs[i][j] = inv * bincoeff(d,j) * bincoeff(t,i-j); end end % } % for i=ph2+1:ph-1 % for (i = ph2+1; i <= ph-1; i++) { mpi = min(d,i); % mpi = min(d, i); for j=max(0,i-t):mpi % for (j = max(0,i-t); j <= mpi; j++) bezalfs(j+1,i+1) = bezalfs(d-j+1,ph-i+1); % bezalfs[i][j] = bezalfs[ph-i][d-j]; end end % } % mh = ph; % mh = ph; kind = ph+1; % kind = ph+1; r = -1; % r = -1; a = d; % a = d; b = d+1; % b = d+1; cind = 1; % cind = 1; ua = k(1); % ua = k[0]; % for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) ic(ii+1,1) = c(ii+1,1); % ictrl[0][ii] = ctrl[0][ii]; end % for i=0:ph % for (i = 0; i <= ph; i++) ik(i+1) = ua; % ik[i] = ua; end % % // initialise first bezier seg for i=0:d % for (i = 0; i <= d; i++) for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) bpts(ii+1,i+1) = c(ii+1,i+1); % bpts[i][ii] = ctrl[i][ii]; end end % % // big loop thru knot vector while b < m % while (b < m) { i = b; % i = b; while b < m && k(b+1) == k(b+2) % while (b < m && k[b] == k[b+1]) b = b + 1; % b++; end % mul = b - i + 1; % mul = b - i + 1; mh = mh + mul + t; % mh += mul + t; ub = k(b+1); % ub = k[b]; oldr = r; % oldr = r; r = d - mul; % r = d - mul; % % // insert knot u(b) r times if oldr > 0 % if (oldr > 0) lbz = floor((oldr+2)/2); % lbz = (oldr+2) / 2; else % else lbz = 1; % lbz = 1; end % if r > 0 % if (r > 0) rbz = ph - floor((r+1)/2); % rbz = ph - (r+1)/2; else % else rbz = ph; % rbz = ph; end % if r > 0 % if (r > 0) { % // insert knot to get bezier segment numer = ub - ua; % numer = ub - ua; for q=d:-1:mul+1 % for (q = d; q > mul; q--) alfs(q-mul) = numer / (k(a+q+1)-ua); % alfs[q-mul-1] = numer / (k[a+q]-ua); end for j=1:r % for (j = 1; j <= r; j++) { save = r - j; % save = r - j; s = mul + j; % s = mul + j; % for q=d:-1:s % for (q = d; q >= s; q--) for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) tmp1 = alfs(q-s+1)*bpts(ii+1,q+1); tmp2 = (1-alfs(q-s+1))*bpts(ii+1,q); bpts(ii+1,q+1) = tmp1 + tmp2; % bpts[q][ii] = alfs[q-s]*bpts[q][ii]+(1.0-alfs[q-s])*bpts[q-1][ii]; end end % for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) Nextbpts(ii+1,save+1) = bpts(ii+1,d+1); % Nextbpts[save][ii] = bpts[d][ii]; end end % } end % } % // end of insert knot % % // degree elevate bezier for i=lbz:ph % for (i = lbz; i <= ph; i++) { for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) ebpts(ii+1,i+1) = 0; % ebpts[i][ii] = 0.0; end mpi = min(d, i); % mpi = min(d, i); for j=max(0,i-t):mpi % for (j = max(0,i-t); j <= mpi; j++) for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) tmp1 = ebpts(ii+1,i+1); tmp2 = bezalfs(j+1,i+1)*bpts(ii+1,j+1); ebpts(ii+1,i+1) = tmp1 + tmp2; % ebpts[i][ii] = ebpts[i][ii] + bezalfs[i][j]*bpts[j][ii]; end end end % } % // end of degree elevating bezier % if oldr > 1 % if (oldr > 1) { % // must remove knot u=k[a] oldr times first = kind - 2; % first = kind - 2; last = kind; % last = kind; den = ub - ua; % den = ub - ua; bet = floor((ub-ik(kind)) / den); % bet = (ub-ik[kind-1]) / den; % % // knot removal loop for tr=1:oldr-1 % for (tr = 1; tr < oldr; tr++) { i = first; % i = first; j = last; % j = last; kj = j - kind + 1; % kj = j - kind + 1; while j-i > tr % while (j - i > tr) { % // loop and compute the new control points % // for one removal step if i < cind % if (i < cind) { alf = (ub-ik(i+1))/(ua-ik(i+1)); % alf = (ub-ik[i])/(ua-ik[i]); for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) tmp1 = alf*ic(ii+1,i+1); tmp2 = (1-alf)*ic(ii+1,i); ic(ii+1,i+1) = tmp1 + tmp2; % ictrl[i][ii] = alf * ictrl[i][ii] + (1.0-alf) * ictrl[i-1][ii]; end end % } if j >= lbz % if (j >= lbz) { if j-tr <= kind-ph+oldr % if (j-tr <= kind-ph+oldr) { gam = (ub-ik(j-tr+1)) / den; % gam = (ub-ik[j-tr]) / den; for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) tmp1 = gam*ebpts(ii+1,kj+1); tmp2 = (1-gam)*ebpts(ii+1,kj+2); ebpts(ii+1,kj+1) = tmp1 + tmp2; % ebpts[kj][ii] = gam*ebpts[kj][ii] + (1.0-gam)*ebpts[kj+1][ii]; end % } else % else { for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) tmp1 = bet*ebpts(ii+1,kj+1); tmp2 = (1-bet)*ebpts(ii+1,kj+2); ebpts(ii+1,kj+1) = tmp1 + tmp2; % ebpts[kj][ii] = bet*ebpts[kj][ii] + (1.0-bet)*ebpts[kj+1][ii]; end end % } end % } i = i + 1; % i++; j = j - 1; % j--; kj = kj - 1; % kj--; end % } % first = first - 1; % first--; last = last + 1; % last++; end % } end % } % // end of removing knot n=k[a] % % // load the knot ua if a ~= d % if (a != d) for i=0:ph-oldr-1 % for (i = 0; i < ph-oldr; i++) { ik(kind+1) = ua; % ik[kind] = ua; kind = kind + 1; % kind++; end end % } % % // load ctrl pts into ic for j=lbz:rbz % for (j = lbz; j <= rbz; j++) { for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) ic(ii+1,cind+1) = ebpts(ii+1,j+1); % ictrl[cind][ii] = ebpts[j][ii]; end cind = cind + 1; % cind++; end % } % if b < m % if (b < m) { % // setup for next pass thru loop for j=0:r-1 % for (j = 0; j < r; j++) for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) bpts(ii+1,j+1) = Nextbpts(ii+1,j+1); % bpts[j][ii] = Nextbpts[j][ii]; end end for j=r:d % for (j = r; j <= d; j++) for ii=0:mc-1 % for (ii = 0; ii < mc; ii++) bpts(ii+1,j+1) = c(ii+1,b-d+j+1); % bpts[j][ii] = ctrl[b-d+j][ii]; end end a = b; % a = b; b = b+1; % b++; ua = ub; % ua = ub; % } else % else % // end knot for i=0:ph % for (i = 0; i <= ph; i++) ik(kind+i+1) = ub; % ik[kind+i] = ub; end end end % } % End big while loop % // end while loop % % *nh = mh - ph - 1; % % freevec2mat(ctrl); % freevec2mat(ictrl); % freematrix(bezalfs); % freematrix(bpts); % freematrix(ebpts); % freematrix(Nextbpts); % mxFree(alfs); % % return(ierr); % } function b = bincoeff(n,k) % Computes the binomial coefficient. % % ( n ) n! % ( ) = -------- % ( k ) k!(n-k)! % % b = bincoeff(n,k) % % Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215. % double bincoeff(int n, int k) % { b = floor(0.5+exp(factln(n)-factln(k)-factln(n-k))); % return floor(0.5+exp(factln(n)-factln(k)-factln(n-k))); % } function f = factln(n) % computes ln(n!) if n <= 1, f = 0; return, end f = gammaln(n+1); %log(factorial(n));</pre>

nurbs_toolbox/bspderiv.m

function [dc,dk] = bspderiv(d,c,k) % % Function Name: % % bspdeval - Evaluate the control points and knot sequence of the derivative % of a univariate B-Spline. % % Calling Sequence: % % [dc,dk] = bspderiv(d,c,k) % % Parameters: % % d : Degree of the B-Spline. % % c : Control Points, matrix of size (dim,nc). % % k : Knot sequence, row vector of size nk. % % dc : Control points of the derivative % % dk : Knot sequence of the derivative % % Description: % % Evaluate the derivative of univariate B-Spline, which is itself a B-Spline. % This function provides an interface to a toolbox 'C' routine. [mc,nc] = size(c); nk = numel(k); % % int bspderiv(int d, double *c, int mc, int nc, double *k, int nk, double *dc, % double *dk) % { % int ierr = 0; % int i, j, tmp; % % // control points % double **ctrl = vec2mat(c,mc,nc); % % // control points of the derivative dc = zeros(mc,nc-1); % double **dctrl = vec2mat(dc,mc,nc-1); % for i=0:nc-2 % for (i = 0; i < nc-1; i++) { tmp = d / (k(i+d+2) - k(i+2)); % tmp = d / (k[i+d+1] - k[i+1]); for j=0:mc-1 % for (j = 0; j < mc; j++) { dc(j+1,i+1) = tmp*(c(j+1,i+2) - c(j+1,i+1)); % dctrl[i][j] = tmp * (ctrl[i+1][j] - ctrl[i][j]); end % } end % } % dk = zeros(1,nk-2); % j = 0; for i=1:nk-2 % for (i = 1; i < nk-1; i++) dk(i) = k(i+1); % dk[j++] = k[i]; end % % freevec2mat(dctrl); % freevec2mat(ctrl); % % return ierr; % }</pre>

nurbs_toolbox/bspeval.m

function p = bspeval(d,c,k,u) % % Function Name: % % bspeval - Evaluate a univariate B-Spline. % % Calling Sequence: % % p = bspeval(d,c,k,u) % % Parameters: % % d : Degree of the B-Spline. % % c : Control Points, matrix of size (dim,nc). % % k : Knot sequence, row vector of size nk. % % u : Parametric evaluation points, row vector of size nu. % % p : Evaluated points, matrix of size (dim,nu) % % Description: % % Evaluate a univariate B-Spline. This function provides an interface to % a toolbox 'C' routine. nu = numel(u); [mc,nc] = size(c); % int bspeval(int d, double *c, int mc, int nc, double *k, int nk, double *u,int nu, double *p){ % int ierr = 0; % int i, s, tmp1, row, col; % double tmp2; % % // Construct the control points % double **ctrl = vec2mat(c,mc,nc); % % // Contruct the evaluated points p = zeros(mc,nu); % double **pnt = vec2mat(p,mc,nu); % % // space for the basis functions N = zeros(d+1,1); % double *N = (double*) mxMalloc((d+1)*sizeof(double)); % % // for each parametric point i for col=1:nu % for (col = 0; col < nu; col++) { % // find the span of u[col] s = findspan(nc-1, d, u(col), k); % s = findspan(nc-1, d, u[col], k); N = basisfun(s,u(col),d,k); % basisfun(s, u[col], d, k, N); % tmp1 = s - d + 1; % tmp1 = s - d; for row=1:mc % for (row = 0; row < mc; row++) { tmp2 = 0; % tmp2 = 0.0; for i=0:d % for (i = 0; i <= d; i++) tmp2 = tmp2 + N(i+1)*c(row,tmp1+i); % tmp2 += N[i] * ctrl[tmp1+i][row]; end % p(row,col) = tmp2; % pnt[col][row] = tmp2; end % } end % } % % mxFree(N); % freevec2mat(pnt); % freevec2mat(ctrl); % % return ierr; % }

nurbs_toolbox/bspkntins.m

function [ic,ik] = bspkntins(d,c,k,u) % % Function Name: % % bspkntins - Insert knots into a univariate B-Spline. % % Calling Sequence: % % [ic,ik] = bspkntins(d,c,k,u) % % Parameters: % % d : Degree of the B-Spline. % % c : Control points, matrix of size (dim,nc). % % k : Knot sequence, row vector of size nk. % % u : Row vector of knots to be inserted, size nu % % ic : Control points of the new B-Spline, of size (dim,nc+nu) % % ik : Knot vector of the new B-Spline, of size (nk+nu) % % Description: % % Insert knots into a univariate B-Spline. This function provides an % interface to a toolbox 'C' routine. [mc,nc] = size(c); nu = numel(u); nk = numel(k); % % int bspkntins(int d, double *c, int mc, int nc, double *k, int nk, % double *u, int nu, double *ic, double *ik) % { % int ierr = 0; % int a, b, r, l, i, j, m, n, s, q, ind; % double alfa; % % double **ctrl = vec2mat(c, mc, nc); ic = zeros(mc,nc+nu); % double **ictrl = vec2mat(ic, mc, nc+nu); ik = zeros(1,nk+nu); % n = size(c,2) - 1; % n = nc - 1; r = length(u) - 1; % r = nu - 1; % m = n + d + 1; % m = n + d + 1; a = findspan(n, d, u(1), k); % a = findspan(n, d, u[0], k); b = findspan(n, d, u(r+1), k); % b = findspan(n, d, u[r], k); b = b+1; % ++b; % for q=0:mc-1 % for (q = 0; q < mc; q++) { for j=0:a-d, ic(q+1,j+1) = c(q+1,j+1); end % for (j = 0; j <= a-d; j++) ictrl[j][q] = ctrl[j][q]; for j=b-1:n, ic(q+1,j+r+2) = c(q+1,j+1); end % for (j = b-1; j <= n; j++) ictrl[j+r+1][q] = ctrl[j][q]; end % } for j=0:a, ik(j+1) = k(j+1); end % for (j = 0; j <= a; j++) ik[j] = k[j]; for j=b+d:m, ik(j+r+2) = k(j+1); end % for (j = b+d; j <= m; j++) ik[j+r+1] = k[j]; % i = b + d - 1; % i = b + d - 1; s = b + d + r; % s = b + d + r; for j=r:-1:0 % for (j = r; j >= 0; j--) { while u(j+1) <= k(i+1) && i > a % while (u[j] <= k[i] && i > a) { for q=0:mc-1 % for (q = 0; q < mc; q++) ic(q+1,s-d) = c(q+1,i-d); % ictrl[s-d-1][q] = ctrl[i-d-1][q]; end ik(s+1) = k(i+1); % ik[s] = k[i]; s = s - 1; % --s; i = i - 1; % --i; end % } for q=0:mc-1 % for (q = 0; q < mc; q++) ic(q+1,s-d) = ic(q+1,s-d+1); % ictrl[s-d-1][q] = ictrl[s-d][q]; end for l=1:d % for (l = 1; l <= d; l++) { ind = s - d + l; % ind = s - d + l; alfa = ik(s+l+1) - u(j+1); % alfa = ik[s+l] - u[j]; if abs(alfa) == 0 % if (fabs(alfa) == 0.0) for q=0:mc-1 % for (q = 0; q < mc; q++) ic(q+1,ind) = ic(q+1,ind+1); % ictrl[ind-1][q] = ictrl[ind][q]; end else % else { alfa = alfa/(ik(s+l+1) - k(i-d+l+1)); % alfa /= (ik[s+l] - k[i-d+l]); for q=0:mc-1 % for (q = 0; q < mc; q++) tmp = (1-alfa)*ic(q+1,ind+1); ic(q+1,ind) = alfa*ic(q+1,ind) + tmp; % ictrl[ind-1][q] = alfa*ictrl[ind-1][q]+(1.0-alfa)*ictrl[ind][q]; end end % } end % } % ik(s+1) = u(j+1); % ik[s] = u[j]; s = s - 1; % --s; end % } % % freevec2mat(ctrl); % freevec2mat(ictrl); % % return ierr; % }

nurbs_toolbox/Contents.m

% NURBS Toolbox. % Version 1.0 % % demos - NURBS demonstrations % % nrbmak - Construct a NURBS from control points and knots. % nrbtform - Applying scaling, translation or rotation operators. % nrbkntins - Knot insertion/refinement. % nrbdegelev - Degree elevation. % nrbderiv - NURBS representation of the derivative. % nrbdeval - Evaluation of the NURBS derivative. % nrbkntmult - Find the multiplilicity of a knot vector. % nrbreverse - Reverse evaluation direction of the NURBS. % nrbtransp - Swap U and V for NURBS surface. % nrbline - Construct a straight line. % nrbcirc - Construct a circular arc. % nrbrect - Construct a rectangle. % nrb4surf - Surface defined by 4 corner points. % nrbeval - Evalution of NURBS curve or surface. % nrbextrude - Extrude a NURBS curve along a vector. % nrbrevolve - Construct surface by revolving a profile. % nrbruled - Ruled surface between twp NURBS curves. % nrbcoons - Construct Coons bilinearly blended surface patch. % nrbplot - Plot NURBS curve or surface. % % bspeval - Evaluate a univariate B-Spline. % bspderiv - B-Spline representation of the derivative % bspkntins - Insert a knot or knots into a univariate B-Spline. % bspdegelev - Degree elevation of a univariate B-Spline. % % vecnorm - Normalise the vectors. % vecmag - Magnitaude of the vectors. % vecmag2 - Squared Magnitude of the vectors. % vecangle - Alternative to atan2 (0 <= angle <= 2*pi) % vecdot - Dot product of two vectors. % veccross - Cross product of two vectors. % vecrotx - Rotation matrix around the x-axis. % vecroty - Rotation matrix around the y-axis. % vecrotz - Rotation matrix around the z-axis. % vecscale - Scaling matrix. % vectrans - Translation matrix. % % deg2rad - Convert degrees to radians. % rad2deg - Convert radians to degrees.

nurbs_toolbox/deg2rad.m

function rad = deg2rad(deg) % % Function Name: % % deg2rad - Convert degrees to radians. % % Calling Sequence: % % rad = deg2rad(deg); % % Parameters: % % deg : Angle in degrees. % % rad : Angle in radians. % % Description: % % Convenient utility function for converting degrees to radians, which are % often the required angular units for functions in the NURBS toolbox. % % Examples: % % // Convert 35 degrees to radians % rad = deg2rad(35); % D.M. Spink % Copyright (c) 2000. rad = pi*deg/180.0;

nurbs_toolbox/demo4surf.m

function demo4surf % Demonstration of a bilinear surface. % % D.M. Spink % Copyright (c) 2000 srf = nrb4surf([0.0 0.0 0.5],[1.0 0.0 -0.5],[0.0 1.0 -0.5],[1.0 1.0 0.5]); nrbplot(srf,[10,10]); title('Construction of a bilinear surface.');

nurbs_toolbox/democirc.m

function democirc % Demonstration of a circle and arcs in the x-y plane. % % D.M. Spink % Copyright (c) 2000 for r = 1:9 crv = nrbcirc(r,[],deg2rad(45),deg2rad(315)); nrbplot(crv,50); hold on; end hold off; axis equal; title('NURBS construction of a 2D circle and arc.');

nurbs_toolbox/democoons.m

% Construction of a bilinearly blended Coons surface % % D.M. Spink % Copyright (c) 2000. % Boundary curve 1 pnts = [ 0.0 3.0 4.5 6.5 8.0 10.0; 0.0 0.0 0.0 0.0 0.0 0.0; 2.0 2.0 7.0 4.0 7.0 9.0]; crv1 = nrbmak(pnts, [0 0 0 1/3 0.5 2/3 1 1 1]); % Boundary curve 2 pnts= [ 0.0 3.0 5.0 8.0 10.0; 10.0 10.0 10.0 10.0 10.0; 3.0 5.0 8.0 6.0 10.0]; crv2 = nrbmak(pnts, [0 0 0 1/3 2/3 1 1 1]); % Boundary curve 3 pnts= [ 0.0 0.0 0.0 0.0; 0.0 3.0 8.0 10.0; 2.0 0.0 5.0 3.0]; crv3 = nrbmak(pnts, [0 0 0 0.5 1 1 1]); % Boundary curve 4 pnts= [ 10.0 10.0 10.0 10.0 10.0; 0.0 3.0 5.0 8.0 10.0; 9.0 7.0 7.0 10.0 10.0]; crv4 = nrbmak(pnts, [0 0 0 0.25 0.75 1 1 1]); srf = nrbcoons(crv1, crv2, crv3, crv4); % Draw the surface nrbplot(srf,[20 20]); title('Construction of a bilinearly blended Coons surface.');

nurbs_toolbox/democurve.m

function democurve % Shows a simple test curve. % % D.M. Spink % Copyright (c) 2000 crv = nrbtestcrv; % plot the control points plot(crv.coefs(1,:),crv.coefs(2,:),'ro'); title('Arbitrary Test 2D Curve.'); hold on; plot(crv.coefs(1,:),crv.coefs(2,:),'r--'); % plot the nurbs curve nrbplot(crv,48); hold off; crv.knots(3)=0.1; figure % plot the control points plot(crv.coefs(1,:),crv.coefs(2,:),'ro'); title('Arbitrary Test 2D Curve.'); hold on; plot(crv.coefs(1,:),crv.coefs(2,:),'r--'); % plot the nurbs curve nrbplot(crv,48); hold off;

nurbs_toolbox/democylind.m

function democylind % Demonstration of the construction of a cylinder. % % D.M. Spink % Copyright (c) 2000 srf = nrbcylind(3,1,[],deg2rad(270),deg2rad(180)); nrbplot(srf,[20,20]); axis equal; title('Cylinderical section by extrusion of a circular arc.');

nurbs_toolbox/demodegelev.m

% Demonstration of the degree elevation algorithm. % crv = nrbtestcrv; % plot the control points plot(crv.coefs(1,:),crv.coefs(2,:),'bo') title('Degree elevation of test curve by 1'); hold on; plot(crv.coefs(1,:),crv.coefs(2,:),'b--'); % draw nurbs curve nrbplot(crv,48); % degree elevate the curve by 1 icrv = nrbdegelev(crv, 1); % insert new knots and plot new control points plot(icrv.coefs(1,:),icrv.coefs(2,:),'ro') plot(icrv.coefs(1,:),icrv.coefs(2,:),'r--'); hold off;

nurbs_toolbox/demodercrv.m

function demodercrv % Demonstrates the construction of a general % curve and determine of the derivative. % % D.M. Spink % Copyright (c) 2000 % make and draw nurbs test curve crv = nrbtestcrv; nrbplot(crv,48); title('First derivatives along a test curve.'); npts = 9; tt = linspace(0.0,1.0,npts); dcrv = nrbderiv(crv); % first derivative [p1, dp] = nrbdeval(crv,dcrv,tt); p2 = vecnorm(dp); hold on; plot(p1(1,:),p1(2,:),'ro'); h = quiver(p1(1,:),p1(2,:),p2(1,:),p2(2,:),0); set(h,'Color','black'); hold off;

nurbs_toolbox/demodersrf.m

function demodersrf % Demonstrates the construction of a general % surface derivatives. % % D.M. Spink % Copyright (c) 2000 % make and draw a test surface srf = nrbtestsrf; p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)}); h = surf(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:))); set(h,'FaceColor','blue','EdgeColor','blue'); title('First derivatives over a test surface.'); npts = 5; tt = linspace(0.0,1.0,npts); dsrf = nrbderiv(srf); [p1, dp] = nrbdeval(srf, dsrf, {tt, tt}); up2 = vecnorm(dp{1}); vp2 = vecnorm(dp{2}); hold on; plot3(p1(1,:),p1(2,:),p1(3,:),'ro'); h1 = quiver3(p1(1,:),p1(2,:),p1(3,:),up2(1,:),up2(2,:),up2(3,:)); h2 = quiver3(p1(1,:),p1(2,:),p1(3,:),vp2(1,:),vp2(2,:),vp2(3,:)); set(h1,'Color','black'); set(h2,'Color','black'); hold off;

nurbs_toolbox/demoellip.m

function demoellip % Demonstration of a unit circle transformed to a inclined ellipse % by first scaling, then rotating and finally translating. % % D.M. Spink % Copyright (c) 2000 xx = vectrans([2.0 1.0])*vecroty(pi/8)*vecrotx(pi/4)*vecscale([1.0 2.0]); c0 = nrbtform(nrbcirc, xx); nrbplot(c0,50); title('Construction of an ellipse by transforming a unit circle.'); grid on;

nurbs_toolbox/demogeom.m

function demogeom % Demonstration of how to construct a 2D geometric % shape from a piece-wise set of NURBSs % % D.M. Spink % Copyright (c) 2000. coefs = [0.0 7.5 15.0 25.0 35.0 30.0 27.5 30.0; 0.0 2.5 0.0 -5.0 5.0 15.0 22.5 30.0]; knots = [0.0 0.0 0.0 1/6 1/3 1/2 2/3 5/6 1.0 1.0 1.0]; % Geometric boundary curve geom = [ nrbmak(coefs,knots) nrbline([30.0 30.0],[20.0 30.0]) nrbline([20.0 30.0],[20.0 20.0]) nrbcirc(10.0,[10.0 20.0],1.5*pi,0.0) nrbline([10.0 10.0],[0.0 10.0]) nrbline([0.0 10.0],[0.0 0.0]) nrbcirc(5.0,[22.5 7.5]) ]; ng = length(geom); for i = 1:ng nrbplot(geom(i),50); hold on; end hold off; axis equal; title('2D Geometry formed by a series of NURBS curves');

nurbs_toolbox/demohelix.m

function demohelix % Demonstration of a 3D helical curve % % D.M. Spink % Copyright (c) 2000 coefs =[ 6.0 0.0 6.0 1; -5.5 0.5 5.5 1; -5.0 1.0 -5.0 1; 4.5 1.5 -4.5 1; 4.0 2.0 4.0 1; -3.5 2.5 3.5 1; -3.0 3.0 -3.0 1; 2.5 3.5 -2.5 1; 2.0 4.0 2.0 1; -1.5 4.5 1.5 1; -1.0 5.0 -1.0 1; 0.5 5.5 -0.5 1; 0.0 6.0 0.0 1]'; knots = [0 0 0 0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1 1 1 1]; crv = nrbmak(coefs,knots); nrbplot(crv,100); title('3D helical curve.'); grid on;

nurbs_toolbox/demokntins.m

function demokntins % Demonstration of the knot insertion algorithm. % % D.M. Spink % Copyright (c) 2000. crv = nrbtestcrv; % plot the control points plot(crv.coefs(1,:),crv.coefs(2,:),'bo') title('Knot insertion along test curve.'); hold on; plot(crv.coefs(1,:),crv.coefs(2,:),'b--'); % draw nurbs curve nrbplot(crv,48); % insert new knots and plot new control points icrv = nrbkntins(crv,[0.125 0.375 0.625 0.875] ); plot(icrv.coefs(1,:),icrv.coefs(2,:),'ro') plot(icrv.coefs(1,:),icrv.coefs(2,:),'r--'); hold off;

nurbs_toolbox/demoline.m

nurbs_toolbox/demorect.m

function demorect % Demonstrate of rectangluar curve % % D.M. Spink % Copyright (c) 2000 crv = nrbtform(nrbrect(2,1), vecrotz(deg2rad(35))); nrbplot(crv,4); title('Construction and rotation of a rectangle.'); axis equal;

nurbs_toolbox/demorevolve.m

function demorevolve % Demonstration of surface construction by revolving a % profile curve. % D.M. Spink % Copyright (c) 2000 % Construct a test profile in the x-z plane pnts = [3.0 5.5 5.5 1.5 1.5 4.0 4.5; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.5 1.5 4.5 3.0 7.5 6.0 8.5]; crv = nrbmak(pnts,[0 0 0 1/4 1/2 3/4 3/4 1 1 1]); % rotate and vectrans by some arbitrary amounts. xx = vecrotz(deg2rad(25))*vecroty(deg2rad(15))*vecrotx(deg2rad(20)); nrb = nrbtform(crv,vectrans([5 5])*xx); % define axes of rotation pnt = [5 5 0]'; vec = xx*[0 0 1 1]'; srf = nrbrevolve(nrb,pnt,vec(1:3)); % make and draw nurbs curve p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)}); surfl(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:))); title('Construct of a 3D surface by revolution of a curve.'); shading interp; colormap(copper); axis equal;

nurbs_toolbox/demoruled.m

% Demonstration of ruled surface construction. % % D.M. Spink % Copyright (c) 2000. % clear all recorded graphics and the current window title('Ruled surface construction from two NURBS curves.'); crv1 = nrbtestcrv; crv2 = nrbtform(nrbcirc(4,[4.5;0],pi,0.0),vectrans([0.0 4.0 -4.0])); srf = nrbruled(crv1,crv2); nrbplot(srf,[40 20]); title('Ruled surface construction from two NURBS curves.');

nurbs_toolbox/demos.m

function tbxStruct=Demos % Demos Demo List infomation for NURBS Toolbox % D.M. Spink % Copyright (c) 2000 if nargout == 0; demo toolbox; return; end tbxStruct.Name='nurbs'; tbxStruct.Type='Toolbox'; tbxStruct.Help= ... {' The NURBS Toolbox provides commands for the construction' ' and use of Non-Uniform Rational B-Splines (NURBS).' ' '}; tbxStruct.DemoList={ '3D Line','demoline', 'Rectangle','demorect', 'Circle and Arc','democirc', 'Helix','demohelix', 'Ellipse', 'demoellip', 'Test Curve','democurve', 'Test Surface','demosurf', 'Bilinear Surface','demo4surf', 'Knot Insertion','demokntins', 'Degree Elevation','demodegelev', 'Curve derivatives','demodercrv', 'Surface derivatives','demodersrf' 'Cylinderical arc','democylind', 'Ruled Surface','demoruled', 'Coons Surface','democoons', 'Test Extrusion','demoextrude', 'Test Revolution - Cup','demorevolve', 'Test Revolution - Ball & Torus','demotorus', '2D Geomtery','demogeom'};

nurbs_toolbox/demotorus.m

function demotorus % A second demonstration of surface construction % by revolution. % % D.M. Spink % Copyright (c) 2000 sphere = nrbrevolve(nrbcirc(1,[],0.0,pi),[0.0 0.0 0.0],[1.0 0.0 0.0]); nrbplot(sphere,[20 20],'light','on'); title('Ball and torus - surface construction by revolution'); hold on; torus = nrbrevolve(nrbcirc(0.2,[0.9 1.0]),[0.0 0.0 0.0],[1.0 0.0 0.0]); nrbplot(torus,[20 10],'light','on'); nrbplot(nrbtform(torus,vectrans([-1.8])),[20 10],'light','on'); hold off;

nurbs_toolbox/findspan.m

function s = findspan(n,p,u,U) % FINDSPAN Find the span of a B-Spline knot vector at a parametric point % ------------------------------------------------------------------------- % ADAPTATION of FINDSPAN from C % ------------------------------------------------------------------------- % % Calling Sequence: % % s = findspan(n,p,u,U) % % INPUT: % % n - number of control points - 1 % p - spline degree % u - parametric point % U - knot sequence % % RETURN: % % s - knot span % % Algorithm A2.1 from 'The NURBS BOOK' pg68 % int findspan(int n, int p, double u, double *U) { % int low, high, mid; % // special case if (u==U(n+2)), s=n; return, end % if (u == U[n+1]) return(n); % % // do binary search low = p; % low = p; high = n + 1; % high = n + 1; mid = floor((low + high) / 2); % mid = (low + high) / 2; while (u < U(mid+1) || u >= U(mid+2)) % while (u < U[mid] || u >= U[mid+1]) { if (u < U(mid+1)) % if (u < U[mid]) high = mid; % high = mid; else % else low = mid; % low = mid; end mid = floor((low + high) / 2); % mid = (low + high) / 2; end % } % s = mid; % return(mid); % }

nurbs_toolbox/gnurbs.m

function out = gnurbs(mode,res) % GNURBS Interactive manipulation of NURBS surface/curve % % For use with the NURBS toolbox, GNURBS allows the user to have an % intuitive manipulation of a NURBS surface/curve via control points. % % GNURBS with no inputs will generate a test surface for manipulation. % GNURBS(NRB), where NRB is a structure defining a NURBS surface/curve % (as defined in the NURBS toolbox). % GNURBS(NRB,RES), where RES is a 2 element vector defines the % resolution of surface/curve rendering in the u and v directions, % ie RES = [URES,VRES]. The default is [20,20] if no argument is % specified. % % <a href="http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=312&objectType=file">NURBS toolbox</a> can be downloaded for free at % http://www.mathworks.com/matlabcentral/fileexchange/ % % Example usage: % % Generate a NURBS surface % gnurbs if ~nargin mode = 'start'; nrb = nrbtestsrf; % Default Surface res = [20 20]; % Default Resolution elseif isstruct(mode) nrb=mode; mode='start'; if nargin<2 res = [20 20]; % Default Resolution end elseif isempty(mode) nrb=nrbtestsrf; mode='start'; if nargin<2 res = [20 20]; % Default Resolution end elseif ischar(mode) fig = gcf; ud = get(gcbo,'userdata'); if strcmp(get(gcbo,'type'),'axes') hnet = get(gcbo,'userdata'); ud = get(hnet,'userdata'); end % disp(mode) else error('GNURBS: Wrong type of input') end switch mode case 'start' % Build surface from NURBS structure [x,y,z] = nrb2xyz(nrb,res); % Get control net [xn,yn,zn] = nrb2net(nrb); % Plot fig = gcf; ax = gca; axis equal n = length(nrb.order); if n > 1 hsurf = surface(x,y,z,'parent',ax,'facecolor','interp',... 'facealpha',.5,'edgealpha',.2); hnet = surface(xn,yn,zn,'parent',ax,'facecolor','none',... 'edgecolor','k','marker','o','linestyle',':',... 'markeredgecolor','k','markerfacecolor','r',... 'markersize',5); else hsurf = line(x,y,z,'parent',ax,'linestyle','-','color','k'); hnet = line(xn,yn,zn,'parent',ax,... 'color','k','marker','o','linestyle',':',... 'markeredgecolor','k','markerfacecolor','r',... 'markersize',5); end hcur = line(0,0,0,'marker','o','markerfacecolor','g',... 'markeredgecolor','k','hittest','off','visible','off',... 'markersize',5); % Store handles and data in figure userdata ud.ax = ax; ud.hsurf = hsurf; ud.hnet = hnet; ud.hcur = hcur; ud.nrb = nrb; ud.res = res; set(hnet,'UserData',ud); % Set Callbacks set(hsurf,'ButtonDownFcn','','hittest','off'); set([hnet gca], 'ButtonDownFcn','gnurbs(''click'')'); set(gca,'userdata',hnet); % Outputs if nargout out = hsurf; end case 'click' [dist,ind] = getClosestPtToLine(ud.hnet); % Don't do anything if too far from closest point if dist > selectThreshold return end ud.ind = ind; [xn,yn,zn] = nrb2net(ud.nrb); setxyz(ud.hcur,xn(ind),yn(ind),zn(ind)); set(ud.hcur,'visible','on'); storeparams(fig,ud); set(fig,'WindowButtonMotionFcn','gnurbs(''move'')'); set(fig,'WindowButtonUpFcn','gnurbs(''up'')'); % set(fig,'UserData',ud); case 'move' stabVector = get(ud.ax,'CurrentPoint'); [xn,yn,zn] = nrb2net(ud.nrb); ind = ud.ind; ptOnViewPlane = [xn(ind),yn(ind),zn(ind)]; [x,y,z] = vectorPlaneIntersect(stabVector,getViewPlane,ptOnViewPlane); xn(ind) = x; yn(ind) = y; zn(ind) = z; coefs = net2nrb(xn,yn,zn); ud.nrb.coefs = coefs; [X,Y,Z] = nrb2xyz(ud.nrb,ud.res); setxyz(ud.hcur,x,y,z); % Update Highlight point setxyz(ud.hsurf,X,Y,Z); % Update Surface setxyz(ud.hnet,xn,yn,zn); % Update Control Net setxyz(ud.hcur,x,y,z); % Update Hightlight point set(fig,'UserData',ud); set(ud.hnet,'userdata',ud); set(ud.hsurf,'userdata',ud.nrb); case 'up' set(ud.hcur,'visible','off'); resetparams(fig); otherwise error('GNURBS: Wrong input string') end function thresh = selectThreshold axesSize = diff([get(gca,'xlim'); get(gca,'ylim'); get(gca,'zlim')],[],2); thresh = sqrt(sum(axesSize.^2))/50; function N = getViewPlane M = view; % Current View matrix up = [0 0 1]'; % Up vector N = (M(1:3,1:3)\up)'; % Normal vector to viewing plane function [X,Y,Z] = vectorPlaneIntersect(vec,N,pt) P0 = vec(1,:); % 1st point on stabbing vector P1 = vec(2,:); % 2nd point on stabbing vector P2 = pt(:)'; % Point on viewing plane (constrains depth ambiguity) % Solve for point where stabbing vector intersects viewing plane u = dot(N,P2-P0)./dot(N,P1-P0); s = (P1-P0)*u + P0; % Separate components X = s(1); Y = s(2); Z = s(3); function [dist,ind] = getClosestPtToLine(h) % GETCLOSESTPTTOLINE Finds the closest point in a data set to a line % In this case, the data set is defined by a surface and the line is % simply the result of get(gca,'CurrentPoint'). % I = GETCLOSESTPTTOLINE(H) returns the index of the point in H closest to a % mouse click in the figure. H is a handle to a surface object. % Daniel Claxton % 02-Mar-2007 % [email protected] x = get(h,'xdata'); y = get(h,'ydata'); z = get(h,'zdata'); if isempty(z), z = x*0; end n = numel(x); ax = get(h,'parent'); cpt = get(ax,'CurrentPoint'); % Find vector in direction of line u = diff(cpt); % Normalize line vector mag = sqrt(sum(u.^2)); mag = max(eps,mag); u = u/mag; % First Point on line (P1) v = cpt(1,:); d = zeros(n,1); for i=1:n pt = [x(i) y(i) z(i)]; % Make vector between P1 and point w = pt - v; % Dot product gives us the magnitude of projection of w onto line f = u*w'; % Multiply magnitude of projection times unit vector u (direction of line) % Shift from orgin to P1 p = f*u + v; % Find shortest distance from point to point on line d(i) = sqrt(sum((pt-p).^2,2)); end [dist,ind] = min(d); function [x,y,z] = nrb2net(srf) % NRB2NET Get control net from NURBS structure x=squeeze(srf.coefs(1,:,:)); y=squeeze(srf.coefs(2,:,:)); z=squeeze(srf.coefs(3,:,:)); function nrb = net2nrb(x,y,z) % NET2NRB Get NURBS coefficients from control net [m,n]=size(x); x = reshape(x,1,m,n); y = reshape(y,1,m,n); z = reshape(z,1,m,n); one = ones(1,m,n); nrb = cat(1,x,y); nrb = cat(1,nrb,z); nrb = cat(1,nrb,one); function [x,y,z] = nrb2xyz(nrb,res) % NRB2XYZ Evaluate NURBS with a given resolution n = length(nrb.order); if n > 1 nu = res(1); nv = res(2); p = nrbeval(nrb,{linspace(0,1,nu),linspace(0,1,nv)}); else p = nrbeval(nrb,linspace(0,1,res(1))); end x = squeeze(p(1,:,:)); y = squeeze(p(2,:,:)); z = squeeze(p(3,:,:)); function setxyz(h,x,y,z) % SETXYZ Set X,Y,Z Data of graphics handle object set(h,'Xdata',x); set(h,'Ydata',y); set(h,'Zdata',z); if strcmp(get(h,'type'),'surface') set(h,'cdata',z/max(z(:))); end function ud = storeparams(fig,ud) ud.udata = get(fig,'userdata'); ud.down = get(fig,'windowbuttondownfcn'); ud.move = get(fig,'windowbuttonmotionfcn'); ud.up = get(fig,'windowbuttonupfcn'); set(fig,'userdata',ud); function resetparams(fig) ud = get(fig,'userdata'); set(fig,'userdata',ud.udata','windowbuttondownfcn',ud.down,... 'windowbuttonmotionfcn',ud.move,'windowbuttonupfcn',ud.up);

nurbs_toolbox/nrb4surf.m

function srf = nrb4surf(p11,p12,p21,p22) % % Function Name: % % nrb4surf - Constructs a NURBS bilinear surface. % % Calling Sequence: % % srf = nrb4surf(p11,p12,p21,p22) % % Parameters: % % p11 : Cartesian coordinate of the lhs bottom corner point. % % p12 : Cartesian coordinate of the rhs bottom corner point. % % p21 : Cartesian coordinate of the lhs top corner point. % % p22 : Cartesian coordinate of the rhs top corner point. % % srf : NURBS bilinear surface, see nrbmak. % % Description: % % Constructs a bilinear surface defined by four coordinates. % % The position of the corner points % % ^ V direction % | % ---------------- % |p21 p22| % | | % | SRF | % | | % |p11 p12| % -------------------> U direction % % D.M. Spink % Copyright (c) 2000. if nargin < 4 error('Four corner points must be defined'); end coefs = [zeros(3,2,2); ones(1,2,2)]; coefs(1:length(p11),1,1) = p11(:); coefs(1:length(p12),1,2) = p12(:); coefs(1:length(p21),2,1) = p21(:); coefs(1:length(p22),2,2) = p22(:); knots = {[0 0 1 1] [0 0 1 1]}; srf = nrbmak(coefs, knots);

nurbs_toolbox/nrbcirc.m

function curve = nrbcirc(radius,center,sang,eang) % % Function Name: % % nrbcirc - Construct a circular arc. % % Calling Sequence: % % crv = nrbarc() % crv = nrbarc(radius) % crv = nrbarc(radius,center) % crv = nrbarc(radius,center,sang,eang) % % Parameters: % % radius : Radius of the circle, default 1.0 % % center : Center of the circle, default (0,0,0) % % sang : Start angle, default 0 degrees % % eang : End angle 360 degrees % % crv : NURBS curve for a circular arc. % % Description: % % Constructs NURBS data structure for a circular arc in the x-y plane. If % no rhs arguments are supplied a unit circle with center (0.0,0.0) is % constructed. % % Angles are defined as positive in the anti-clockwise direction. % D.M. Spink % Copyright (c) 2000. if nargin < 1 radius = 1; end if nargin < 2 center = []; end if nargin < 4 sang = 0; eang = 2*pi; end sweep = eang - sang; % sweep angle of arc if sweep < 0 sweep = 2*pi + sweep; end if abs(sweep) <= pi/2 narcs = 1; % number of arc segments knots = [0 0 0 1 1 1]; elseif abs(sweep) <= pi narcs = 2; knots = [0 0 0 0.5 0.5 1 1 1]; elseif abs(sweep) <= 3*pi/2 narcs = 3; knots = [0 0 0 1/3 1/3 2/3 2/3 1 1 1]; else narcs = 4; knots = [0 0 0 0.25 0.25 0.5 0.5 0.75 0.75 1 1 1]; end dsweep = sweep/(2*narcs); % arc segment sweep angle/2 % determine middle control point and weight wm = cos(dsweep); x = radius*wm; y = radius*sin(dsweep); xm = x+y*tan(dsweep); % arc segment control points ctrlpt = [ x wm*xm x; % w*x - coordinate -y 0 y; % w*y - coordinate 0 0 0; % w*z - coordinate 1 wm 1]; % w - coordinate % build up complete arc from rotated segments coefs = zeros(4,2*narcs+1); % nurbs control points of arc xx = vecrotz(sang + dsweep); coefs(:,1:3) = xx*ctrlpt; % rotate to start angle xx = vecrotz(2*dsweep); for n = 2:narcs m = 2*n+[0 1]; coefs(:,m) = xx*coefs(:,m-2); end % vectrans arc if necessary if ~isempty(center) xx = vectrans(center); coefs = xx*coefs; end curve = nrbmak(coefs,knots);

nurbs_toolbox/nrbcoons.m

function srf = nrbcoons(u1, u2, v1, v2) % % Function Name: % % nrbcoons - Construction of a Coons patch. % % Calling Sequence: % % srf = nrbcoons(ucrv1, ucrv2, vcrv1, vcrv2) % % Parameters: % % ucrv1 : NURBS curve defining the bottom U direction boundary of % the constructed NURBS surface. % % ucrv2 : NURBS curve defining the top U direction boundary of % the constructed NURBS surface. % % vcrv1 : NURBS curve defining the bottom V direction boundary of % the constructed NURBS surface. % % vcrv2 : NURBS curve defining the top V direction boundary of % the constructed NURBS surface. % % srf : Coons NURBS surface patch. % % Description: % % Construction of a bilinearly blended Coons surface patch from four NURBS % curves that define the boundary. % % The orientation of the four NURBS boundary curves. % % ^ V direction % | % | ucrv2 % ------->-------- % | | % | | % vcrv1 ^ Surface ^ vcrv2 % | | % | | % ------->-----------> U direction % ucrv1 % % % Examples: % % // Define four NURBS curves and construct a Coons surface patch. % pnts = [ 0.0 3.0 4.5 6.5 8.0 10.0; % 0.0 0.0 0.0 0.0 0.0 0.0; % 2.0 2.0 7.0 4.0 7.0 9.0]; % crv1 = nrbmak(pnts, [0 0 0 1/3 0.5 2/3 1 1 1]); % % pnts= [ 0.0 3.0 5.0 8.0 10.0; % 10.0 10.0 10.0 10.0 10.0; % 3.0 5.0 8.0 6.0 10.0]; % crv2 = nrbmak(pnts, [0 0 0 1/3 2/3 1 1 1]); % % pnts= [ 0.0 0.0 0.0 0.0; % 0.0 3.0 8.0 10.0; % 2.0 0.0 5.0 3.0]; % crv3 = nrbmak(pnts, [0 0 0 0.5 1 1 1]); % % pnts= [ 10.0 10.0 10.0 10.0 10.0; % 0.0 3.0 5.0 8.0 10.0; % 9.0 7.0 7.0 10.0 10.0]; % crv4 = nrbmak(pnts, [0 0 0 0.25 0.75 1 1 1]); % % srf = nrbcoons(crv1, crv2, crv3, crv4); % nrbplot(srf,[20 20],220,45); % D.M. Spink % Copyright (c) 2000. if nargin ~= 4 error('Incorrect number of input arguments'); end r1 = nrbruled(u1, u2); r2 = nrbtransp(nrbruled(v1, v2)); t = nrb4surf(u1.coefs(:,1), u1.coefs(:,end), u2.coefs(:,1), u2.coefs(:,end)); % raise all surfaces to a common degree du = max([r1.order(1), r2.order(1), t.order(1)]); dv = max([r1.order(2), r2.order(2), t.order(2)]); r1 = nrbdegelev(r1, [du - r1.order(1), dv - r1.order(2)]); r2 = nrbdegelev(r2, [du - r2.order(1), dv - r2.order(2)]); t = nrbdegelev(t, [du - t.order(1), dv - t.order(2)]); % merge the knot vectors, to obtain a common knot vector % U knots k1 = r1.knots{1}; k2 = r2.knots{1}; k3 = t.knots{1}; k = unique([k1 k2 k3]); n = length(k); kua = []; kub = []; kuc = []; for i = 1:n i1 = length(find(k1 == k(i))); i2 = length(find(k2 == k(i))); i3 = length(find(k3 == k(i))); m = max([i1, i2, i3]); kua = [kua k(i)*ones(1,m-i1)]; kub = [kub k(i)*ones(1,m-i2)]; kuc = [kuc k(i)*ones(1,m-i3)]; end % V knots k1 = r1.knots{2}; k2 = r2.knots{2}; k3 = t.knots{2}; k = unique([k1 k2 k3]); n = length(k); kva = []; kvb = []; kvc = []; for i = 1:n i1 = length(find(k1 == k(i))); i2 = length(find(k2 == k(i))); i3 = length(find(k3 == k(i))); m = max([i1, i2, i3]); kva = [kva k(i)*ones(1,m-i1)]; kvb = [kvb k(i)*ones(1,m-i2)]; kvc = [kvc k(i)*ones(1,m-i3)]; end r1 = nrbkntins(r1, {kua, kva}); r2 = nrbkntins(r2, {kub, kvb}); t = nrbkntins(t, {kuc, kvc}); % combine coefficient to construct Coons surface coefs(1,:,:) = r1.coefs(1,:,:) + r2.coefs(1,:,:) - t.coefs(1,:,:); coefs(2,:,:) = r1.coefs(2,:,:) + r2.coefs(2,:,:) - t.coefs(2,:,:); coefs(3,:,:) = r1.coefs(3,:,:) + r2.coefs(3,:,:) - t.coefs(3,:,:); coefs(4,:,:) = r1.coefs(4,:,:) + r2.coefs(4,:,:) - t.coefs(4,:,:); srf = nrbmak(coefs, r1.knots);

nurbs_toolbox/nrbcylind.m

function surf = nrbcylind(height,radius,center,sang,eang) % % Function Name: % % nrbcylind - Construct a cylinder or cylindrical patch. % % Calling Sequence: % % srf = nrbcylind() % srf = nrbcylind(height) % srf = nrbcylind(height,radius) % srf = nrbcylind(height,radius,center) % srf = nrbcylind(height,radius,center,sang,eang) % % Parameters: % % height : Height of the cylinder along the axis, default 1.0 % % radius : Radius of the cylinder, default 1.0 % % center : Center of the cylinder, default (0,0,0) % % sang : Start angle relative to the origin, default 0. % % eang : End angle relative to the origin, default 2*pi. % % % Description: % % Construct a cylinder or cylindrical patch by extruding a circular arc. % D.M. Spink % Copyright (c) 2000. if nargin < 1 height = 1; end if nargin < 2 radius = 1; end if nargin < 3 center = []; end if nargin < 5 sang = 0; eang = 2*pi; end surf = nrbextrude(nrbcirc(radius,center,sang,eang),[0.0 0.0 height]);

nurbs_toolbox/nrbdegelev.m

function inurbs = nrbdegelev(nurbs, ntimes) % % Function Name: % % nrbdegelev - Elevate the degree of the NURBS curve or surface. % % Calling Sequence: % % ecrv = nrbdegelev(crv,utimes); % esrf = nrbdegelev(srf,{utimes,vtimes}); % % Parameters: % % crv : NURBS curve, see nrbmak. % % srf : NURBS surface, see nrbmak. % % utimes : Increase the degree along U direction utimes. % % vtimes : Increase the degree along V direction vtimes. % % ecrv : new NURBS structure for a curve with degree elevated. % % esrf : new NURBS structure for a surface with degree elevated. % % % Description: % % Degree elevates the NURBS curve or surface. This function uses the % B-Spline function bspdegelev, which interface to an internal 'C' % routine. % % Examples: % % Increase the NURBS surface twice along the V direction. % esrf = nrbdegelev(srf, 0, 1); % % See: % % bspdegelev % D.M. Spink % Copyright (c) 2000. if nargin < 2 error('Input argument must include the NURBS and degree increment.'); end if ~isstruct(nurbs) error('NURBS representation is not structure!'); end if ~strcmp(nurbs.form,'B-NURBS') error('Not a recognised NURBS representation'); end degree = nurbs.order-1; if iscell(nurbs.knots) % NURBS represents a surface [dim,num1,num2] = size(nurbs.coefs); % Degree elevate along the v direction if ntimes(2) == 0 coefs = nurbs.coefs; knots{2} = nurbs.knots{2}; else coefs = reshape(nurbs.coefs,4*num1,num2); [coefs,knots{2}] = bspdegelev(degree(2),coefs,nurbs.knots{2},ntimes(2)); num2 = size(coefs,2); coefs = reshape(coefs,[4 num1 num2]); end % Degree elevate along the u direction if ntimes(1) == 0 knots{1} = nurbs.knots{1}; else coefs = permute(coefs,[1 3 2]); coefs = reshape(coefs,4*num2,num1); [coefs,knots{1}] = bspdegelev(degree(1),coefs,nurbs.knots{1},ntimes(1)); coefs = reshape(coefs,[4 num2 size(coefs,2)]); coefs = permute(coefs,[1 3 2]); end else % NURBS represents a curve if isempty(ntimes) coefs = nurbs.coefs; knots = nurbs.knots; else [coefs,knots] = bspdegelev(degree,nurbs.coefs,nurbs.knots,ntimes); end end % construct new NURBS inurbs = nrbmak(coefs,knots);

nurbs_toolbox/nrbderiv.m

function dnurbs = nrbderiv(nurbs) % % Function Name: % % nrbderiv - Construct the first derivative representation of a % NURBS curve or surface. % % Calling Sequence: % % ders = nrbderiv(nrb); % % Parameters: % % nrb : NURBS data structure, see nrbmak. % % ders : A data structure that represents the first % derivatives of a NURBS curve or surface. % % Description: % % The derivatives of a B-Spline are themselves a B-Spline of lower degree, % giving an efficient means of evaluating multiple derivatives. However, % although the same approach can be applied to NURBS, the situation for % NURBS is a more complex. I have at present restricted the derivatives % to just the first. I don't claim that this implentation is % the best approach, but it will have to do for now. The function returns % a data struture that for a NURBS curve contains the first derivatives of % the B-Spline representation. Remember that a NURBS curve is represent by % a univariate B-Spline using the homogeneous coordinates. % The derivative data structure can be evaluated later with the function % nrbdeval. % % Examples: % % See the function nrbdeval for an example. % % See: % % nrbdeval % D.M. Spink % Copyright (c) 2000. if ~isstruct(nurbs) error('NURBS representation is not structure!'); end if ~strcmp(nurbs.form,'B-NURBS') error('Not a recognised NURBS representation'); end degree = nurbs.order - 1; if iscell(nurbs.knots) % NURBS structure represents a surface num1 = nurbs.number(1); num2 = nurbs.number(2); % taking derivatives along the u direction dknots = nurbs.knots; dcoefs = permute(nurbs.coefs,[1 3 2]); dcoefs = reshape(dcoefs,4*num2,num1); [dcoefs,dknots{1}] = bspderiv(degree(1),dcoefs,nurbs.knots{1}); dcoefs = permute(reshape(dcoefs,[4 num2 size(dcoefs,2)]),[1 3 2]); dnurbs{1} = nrbmak(dcoefs, dknots); % taking derivatives along the v direction dknots = nurbs.knots; dcoefs = reshape(nurbs.coefs,4*num1,num2); [dcoefs,dknots{2}] = bspderiv(degree(2),dcoefs,nurbs.knots{2}); dcoefs = reshape(dcoefs,[4 num1 size(dcoefs,2)]); dnurbs{2} = nrbmak(dcoefs, dknots); else % NURBS structure represents a curve [dcoefs,dknots] = bspderiv(degree,nurbs.coefs,nurbs.knots); dnurbs = nrbmak(dcoefs, dknots); end

nurbs_toolbox/nrbdeval.m

function [pnt,jac] = nrbdeval(nurbs, dnurbs, tt) % Evaluation of the derivative NURBS curve or surface. % % [pnt, jac] = nrbdeval(crv, dcrv, tt) % [pnt, jac] = nrbdeval(srf, dsrf, {tu tv}) % % INPUTS: % % crv - original NURBS curve. % % srf - original NUBRS surface % % dcrv - NURBS derivative represention of crv % % dsrf - NURBS derivative represention of surface % % tt - parametric evaluation points % If the nurbs is a surface then tt is a cell % {tu, tv} are the parametric coordinates % % pnt - evaluated points. % jac - evaluated first derivatives (Jacobian). % % Examples: % % // Determine the first derivatives a NURBS curve at 9 points for 0.0 to % // 1.0 % tt = linspace(0.0, 1.0, 9); % dcrv = nrbderiv(crv); % [pnts,jac] = nrbdeval(crv, dcrv, tt); % D.M. Spink % Copyright (c) 2000. if ~isstruct(nurbs) error('NURBS representation is not structure!'); end if ~strcmp(nurbs.form,'B-NURBS') error('Not a recognised NURBS representation'); end [cp,cw] = nrbeval(nurbs, tt); if iscell(nurbs.knots) % NURBS structure represents a surface temp = cw(ones(3,1),:,:); pnt = cp./temp; [cup,cuw] = nrbeval(dnurbs{1}, tt); tempu = cuw(ones(3,1),:,:); jac{1} = (cup-tempu.*pnt)./temp; [cvp,cvw] = nrbeval(dnurbs{2}, tt); tempv = cvw(ones(3,1),:,:); jac{2} = (cvp-tempv.*pnt)./temp; else % NURBS is a curve temp = cw(ones(3,1),:); pnt = cp./temp; % first derivative [cup,cuw] = nrbeval(dnurbs,tt); temp1 = cuw(ones(3,1),:); jac = (cup-temp1.*pnt)./temp; end

nurbs_toolbox/nrbeval.m

function [p,w] = nrbeval(nurbs,tt) % % Function Name: % % nrbeval - Evaluate a NURBS at parameteric points % % Calling Sequence: % % [p,w] = nrbeval(crv,ut) % [p,w] = nrbeval(srf,{ut,vt}) % % Parameters: % % crv : NURBS curve, see nrbmak. % % srf : NURBS surface, see nrbmak. % % ut : Parametric evaluation points along U direction. % % vt : Parametric evaluation points along V direction. % % p : Evaluated points on the NURBS curve or surface as cartesian % coordinates (x,y,z). If w is included on the lhs argument list % the points are returned as homogeneous coordinates (wx,wy,wz). % % w : Weights of the homogeneous coordinates of the evaluated % points. Note inclusion of this argument changes the type % of coordinates returned in p (see above). % % Description: % % Evaluation of NURBS curves or surfaces at parametric points along the % U and V directions. Either homogeneous coordinates are returned if the % weights are requested in the lhs arguments, or as cartesian coordinates. % This function utilises the 'C' interface bspeval. % % Examples: % % Evaluate the NURBS circle at twenty points from 0.0 to 1.0 % % nrb = nrbcirc; % ut = linspace(0.0,1.0,20); % p = nrbeval(nrb,ut); % % See: % % bspeval % % D.M. Spink % Copyright (c) 2000. if nargin < 2 error('Not enough input arguments'); end foption = 1; % output format 3D cartesian coordinates if nargout == 2 foption = 0; % output format 4D homogenous coordinates end if ~isstruct(nurbs) error('NURBS representation is not structure!'); end if ~strcmp(nurbs.form,'B-NURBS') error('Not a recognised NURBS representation'); end if iscell(nurbs.knots) % NURBS structure represents a surface num1 = nurbs.number(1); num2 = nurbs.number(2); degree = nurbs.order-1; if iscell(tt) % Evaluate over a [u,v] grid % tt{1} represents the u direction % tt{2} represents the v direction nt1 = length(tt{1}); nt2 = length(tt{2}); % Evaluate along the v direction val = reshape(nurbs.coefs,4*num1,num2); val = bspeval(degree(2),val,nurbs.knots{2},tt{2}); val = reshape(val,[4 num1 nt2]); % Evaluate along the u direction val = permute(val,[1 3 2]); val = reshape(val,4*nt2,num1); val = bspeval(degree(1),val,nurbs.knots{1},tt{1}); val = reshape(val,[4 nt2 nt1]); val = permute(val,[1 3 2]); w = val(4,:,:); p = val(1:3,:,:); if foption p = p./repmat(w,[3 1 1]); end else % Evaluate at scattered points % tt(1,:) represents the u direction % tt(2,:) represents the v direction nt = size(tt,2); val = reshape(nurbs.coefs,4*num1,num2); val = bspeval(degree(2),val,nurbs.knots{2},tt(2,:)); val = reshape(val,[4 num1 nt]); % evaluate along the u direction pnts = zeros(4,nt); for v = 1:nt coefs = squeeze(val(:,:,v)); pnts(:,v) = bspeval(degree(1),coefs,nurbs.knots{1},tt(1,v)); end w = pnts(4,:); p = pnts(1:3,:); if foption p = p./w; end end else % NURBS structure represents a curve % tt represent a vector of parametric points in the u direction val = bspeval(nurbs.order-1,nurbs.coefs,nurbs.knots,tt); w = val(4,:); p = val(1:3,:); if foption p = p./repmat(w,3,1); end end

nurbs_toolbox/nrbextrude.m

function srf = nrbextrude(curve,vector) % % Function Name: % % nrbextrude - Construct a NURBS surface by extruding a NURBS curve. % % Calling Sequence: % % srf = nrbextrude(crv,vec); % % Parameters: % % crv : NURBS curve to extrude, see nrbmak. % % vec : Vector along which the curve is extruded. % % srf : NURBS surface constructed. % % Description: % % Constructs a NURBS surface by extruding a NURBS curve along a defined % vector. The NURBS curve forms the U direction of the surface edge, and % extruded along the vector in the V direction. Note NURBS surfaces cannot % be extruded. % % Examples: % % Form a hollow cylinder by extruding a circle along the z-axis. % % srf = nrbextrude(nrbcirc, [0,0,1]); % D.M. Spink % Copyright (c) 2000. if iscell(curve.knots) error('Nurb surfaces cannot be extruded!'); end if nargin < 2 error('Error too few input arguments!'); end if nargin == 3 dz = 0.0; end coefs = cat(3,curve.coefs,vectrans(vector)*curve.coefs); srf = nrbmak(coefs,{curve.knots [0 0 1 1]});

nurbs_toolbox/nrbkntins.m

function inurbs = nrbkntins(nurbs,iknots) % % Function Name: % % nrbkntins - Insert a single or multiple knots into a NURBS curve or % surface. % % Calling Sequence: % % icrv = nrbkntins(crv,iuknots); % isrf = nrbkntins(srf,{iuknots ivknots}); % % Parameters: % % crv : NURBS curve, see nrbmak. % % srf : NURBS surface, see nrbmak. % % iuknots : Knots to be inserted along U direction. % % ivknots : Knots to be inserted along V direction. % % icrv : new NURBS structure for a curve with knots inserted. % % isrf : new NURBS structure for a surface with knots inserted. % % Description: % % Inserts knots into the NURBS data structure, these can be knots at % new positions or at the location of existing knots to increase the % multiplicity. Note that the knot multiplicity cannot be increased % beyond the order of the spline. Knots along the V direction can only % inserted into NURBS surfaces, not curve that are always defined along % the U direction. This function use the B-Spline function bspkntins, % which interfaces to an internal 'C' routine. % % Examples: % % Insert two knots into a curve, one at 0.3 and another % twice at 0.4 % % icrv = nrbkntins(crv, [0.3 0.4 0.4]) % % Insert into a surface two knots as (1) into the U knot % sequence and one knot into the V knot sequence at 0.5. % % isrf = nrbkntins(srf, {[0.3 0.4 0.4] [0.5]}) % % Also see: % % bspkntins % % Note: % % No knot multiplicity will be increased beyond the order of the spline. % D.M. Spink % Copyright (c) 2000. if nargin < 2 error('Input argument must include the NURBS and knots to be inserted'); end if ~isstruct(nurbs) error('NURBS representation is not structure!'); end if ~strcmp(nurbs.form,'B-NURBS') error('Not a recognised NURBS representation'); end degree = nurbs.order-1; if iscell(nurbs.knots) % NURBS represents a surface num1 = nurbs.number(1); num2 = nurbs.number(2); % Insert knots along the v direction if isempty(iknots{2}) coefs = nurbs.coefs; knots{2} = nurbs.knots{2}; else coefs = reshape(nurbs.coefs,4*num1,num2); [coefs,knots{2}] = bspkntins(degree(2),coefs,nurbs.knots{2},iknots{2}); num2 = size(coefs,2); coefs = reshape(coefs,[4 num1 num2]); end % Insert knots along the u direction if isempty(iknots{1}) knots{1} = nurbs.knots{1}; else coefs = permute(coefs,[1 3 2]); coefs = reshape(coefs,4*num2,num1); [coefs,knots{1}] = bspkntins(degree(1),coefs,nurbs.knots{1},iknots{1}); coefs = reshape(coefs,[4 num2 size(coefs,2)]); coefs = permute(coefs,[1 3 2]); end else % NURBS represents a curve if isempty(iknots) coefs = nurbs.coefs; knots = nurbs.knots; else [coefs,knots] = bspkntins(degree,nurbs.coefs,nurbs.knots,iknots); end end % construct new NURBS inurbs = nrbmak(coefs,knots);

nurbs_toolbox/nrbline.m

function curve = nrbline(p1,p2) % % Function Name: % % nrbline - Construct a straight line. % % Calling Sequence: % % crv = nrbline() % crv = nrbline(p1,p2) % % Parameters: % % p1 : 2D or 3D cartesian coordinate of the start point. % % p2 : 2D or 3D cartesian coordinate of the end point. % % crv : NURBS curve for a straight line. % % Description: % % Constructs NURBS data structure for a straight line. If no rhs % coordinates are included the function returns a unit straight % line along the x-axis. % D.M. Spink % Copyright (c) 2000. coefs = [zeros(3,2); ones(1,2)]; if nargin < 2 coefs(1,2) = 1.0; else coefs(1:length(p1),1) = p1(:); coefs(1:length(p2),2) = p2(:); end curve = nrbmak(coefs, [0 0 1 1]);

nurbs_toolbox/nrbmak.m

function nurbs = nrbmak(coefs,knots) % % Function Name: % % nrbmak - Construct the NURBS structure given the control points % and the knots. % % Calling Sequence: % % nurbs = nrbmak(cntrl,knots); % % Parameters: % % cntrl : Control points, these can be either Cartesian or % homogeneous coordinates. % % For a curve the control points are represented by a % matrix of size (dim,nu) and for a surface a multidimensional % array of size (dim,nu,nv). Where nu is number of points along % the parametric U direction, and nv the number of points % along the V direction. Dim is the dimension valid options % are % 2 .... (x,y) 2D Cartesian coordinates % 3 .... (x,y,z) 3D Cartesian coordinates % 4 .... (wx,wy,wz,w) 4D homogeneous coordinates % % knots : Non-decreasing knot sequence spanning the interval % [0.0,1.0]. It's assumed that the curves and surfaces % are clamped to the start and end control points by knot % multiplicities equal to the spline order. % For curve knots form a vector and for a surface the knot % are stored by two vectors for U and V in a cell structure. % {uknots vknots} % % nurbs : Data structure for representing a NURBS curve. % % NURBS Structure: % % Both curves and surfaces are represented by a structure that is % compatible with the Spline Toolbox from Mathworks % % nurbs.form .... Type name 'B-NURBS' % nurbs.dim .... Dimension of the control points % nurbs.number .... Number of Control points % nurbs.coefs .... Control Points % nurbs.order .... Order of the spline % nurbs.knots .... Knot sequence % % Note: the control points are always converted and stored within the % NURBS structure as 4D homogeneous coordinates. A curve is always stored % along the U direction, and the vknots element is an empty matrix. For % a surface the spline degree is a vector [du,dv] containing the degree % along the U and V directions respectively. % % Description: % % This function is used as a convenient means of constructing the NURBS % data structure. Many of the other functions in the toolbox rely on the % NURBS structure been correctly defined as shown above. The nrbmak not % only constructs the proper structure, but also checks for consistency. % The user is still free to build his own structure, in fact a few % functions in the toolbox do this for convenience. % % Examples: % % Construct a 2D line from (0.0,0.0) to (1.5,3.0). % For a straight line a spline of order 2 is required. % Note that the knot sequence has a multiplicity of 2 at the % start (0.0,0.0) and end (1.0 1.0) in order to clamp the ends. % % line = nrbmak([0.0 1.5; 0.0 3.0],[0.0 0.0 1.0 1.0]); % nrbplot(line, 2); % % Construct a surface in the x-y plane i.e % % ^ (0.0,1.0) ------------ (1.0,1.0) % | | | % | V | | % | | Surface | % | | | % | | | % | (0.0,0.0) ------------ (1.0,0.0) % | % |------------------------------------> % U % % coefs = cat(3,[0 0; 0 1],[1 1; 0 1]); % knots = {[0 0 1 1] [0 0 1 1]} % plane = nrbmak(coefs,knots); % nrbplot(plane, [2 2]); % D.M. Spink % Copyright (c) 2000. nurbs.form = 'B-NURBS'; nurbs.dim = 4; np = size(coefs); dim = np(1); if iscell(knots) % constructing a surface nurbs.number = np(2:3); if (dim < 4) nurbs.coefs = repmat([0.0 0.0 0.0 1.0]',[1 np(2:3)]); nurbs.coefs(1:dim,:,:) = coefs; else nurbs.coefs = coefs; end uorder = size(knots{1},2)-np(2); vorder = size(knots{2},2)-np(3); uknots = sort(knots{1}); vknots = sort(knots{2}); uknots = (uknots-uknots(1))/(uknots(end)-uknots(1)); vknots = (vknots-vknots(1))/(vknots(end)-vknots(1)); nurbs.knots = {uknots vknots}; nurbs.order = [uorder vorder]; else % constructing a curve nurbs.number = np(2); if (dim < 4) nurbs.coefs = repmat([0.0 0.0 0.0 1.0]',[1 np(2)]); nurbs.coefs(1:dim,:) = coefs; else nurbs.coefs = coefs; end nurbs.order = size(knots,2)-np(2); knots = sort(knots); nurbs.knots = (knots-knots(1))/(knots(end)-knots(1)); end

nurbs_toolbox/nrbplot.m

nurbs_toolbox/nrbrect.m

function curve = nrbrect(w,h) % % Function Name: % % nrbrect - Construct NURBS representation of a rectangle. % % Calling Sequence: % % crv = nrbrect() % crv = nrbrect(size) % crv = nrbrect(width, height) % % Parameters: % % size : Size of the square (width = height). % % width : Width of the rectangle (along x-axis). % % height : Height of the rectangle (along y-axis). % % crv : NURBS curve, see nrbmak. % % % Description: % % Construct a rectangle or square in the x-y plane with the bottom % lhs corner at (0,0,0). If no rhs arguments provided the function % constructs a unit square. % D.M. Spink % Copyright (c) 2000. if nargin < 1 w = 1; h = 1; end if nargin < 2 h = w; end coefs = [0 w w w w 0 0 0; 0 0 0 h h h h 0; 0 0 0 0 0 0 0 0; 1 1 1 1 1 1 1 1]; knots = [0 0 0.25 0.25 0.5 0.5 0.75 0.75 1 1]; curve = nrbmak(coefs, knots);

nurbs_toolbox/nrbreverse.m

function rnrb = nrbreverse(nrb) % % Function Name: % % nrbreverse - Reverse the evaluation direction of a NURBS curve or surface. % % Calling Sequence: % % rnrb = nrbreverse(nrb); % % Parameters: % % nrb : NURBS data structure, see nrbmak. % % rnrb : Reversed NURBS. % % Description: % % Utility function to reverse the evaluation direction of a NURBS % curve or surface. % D.M. Spink % Copyright (c) 2000. if nargin ~= 1 error('Incorrect number of input arguments'); end if iscell(nrb.knots) % reverse a NURBS surface coefs = nrb.coefs(:,:,end:-1:1); rnrb = nrbmak(coefs(:,end:-1:1,:), {1.0-fliplr(nrb.knots{1}),... 1.0-fliplr(nrb.knots{2})}); else % reverse a NURBS curve rnrb = nrbmak(fliplr(nrb.coefs), 1.0-fliplr(nrb.knots)); end

nurbs_toolbox/nrbrevolve.m

function surf = nrbrevolve(curve,pnt,vec,theta) % % Function Name: % % nrbrevolve - Construct a NURBS surface by revolving a NURBS curve. % % Calling Sequence: % % srf = nrbrevolve(crv,pnt,vec[,ang]) % % Parameters: % % crv : NURBS curve to revolve, see nrbmak. % % pnt : Coordinate of the point used to define the axis % of rotation. % % vec : Vector defining the direction of the rotation axis. % % ang : Angle to revolve the curve, default 2*pi % % Description: % % Construct a NURBS surface by revolving the profile NURBS curve around % an axis defined by a point and vector. % % Examples: % % Construct a sphere by rotating a semicircle around a x-axis. % % crv = nrbcirc(1.0,[0 0 0],0,pi); % srf = nrbrevolve(crv,[0 0 0],[1 0 0]); % nrbplot(srf,[20 20]); % % NOTE: % % The algorithm: % % 1) vectrans the point to the origin (0,0,0) % 2) rotate the vector into alignment with the z-axis % % for each control point along the curve % % 3) determine the radius and angle of control % point to the z-axis % 4) construct a circular arc in the x-y plane with % this radius and start angle and sweep angle theta % 5) combine the arc and profile, coefs and weights. % % next control point % % 6) rotate and vectrans the surface back into position % by reversing 1 and 2. % % D.M. Spink % Copyright (c) 2000. if nargin < 3 error('Not enough arguments to construct revolved surface'); end if nargin < 4 theta = 2.0*pi; end % Translate curve the center point to the origin if isempty(pnt) pnt = zeros(3,1); end if length(pnt) ~= 3 error('All point and vector coordinates must be 3D'); end % Translate and rotate the curve into alignment with the z-axis T = vectrans(-pnt); angx = vecangle(vec(1),vec(3)); RY = vecroty(-angx); vectmp = RY*[vecnorm(vec(:));1.0]; angy = vecangle(vectmp(2),vectmp(3)); RX = vecrotx(angy); curve = nrbtform(curve,RX*RY*T); % Construct an arc arc = nrbcirc(1.0,[],0.0,theta); % Construct the surface coefs = zeros(4,arc.number,curve.number); angle = vecangle(curve.coefs(2,:),curve.coefs(1,:)); radius = vecmag(curve.coefs(1:2,:)); for i = 1:curve.number coefs(:,:,i) = vecrotz(angle(i))*vectrans([0.0 0.0 curve.coefs(3,i)])*... vecscale([radius(i) radius(i)])*arc.coefs; coefs(4,:,i) = coefs(4,:,i)*curve.coefs(4,i); end surf = nrbmak(coefs,{arc.knots curve.knots}); % Rotate and vectrans the surface back into position T = vectrans(pnt); RX = vecrotx(-angy); RY = vecroty(angx); surf = nrbtform(surf,T*RY*RX);

nurbs_toolbox/nrbruled.m

function srf = nrbruled(crv1, crv2) % % Function Name: % % nrbruled - Constructs a ruled surface between two NURBS curves. % % Calling Sequence: % % srf = nrbruled(crv1, crv2) % % Parameters: % % crv1 : First NURBS curve, see nrbmak. % % crv2 : Second NURBS curve, see nrbmak. % % srf : Ruled NURBS surface. % % Description: % % Constructs a ruled surface between two NURBS curves. The ruled surface is % ruled along the V direction. % % Examples: % % Construct a ruled surface between a semicircle and a straight line. % % cir = nrbcirc(1,[0 0 0],0,pi); % line = nrbline([-1 0.5 1],[1 0.5 1]); % srf = nrbruled(cir,line); % nrbplot(srf,[20 20]); % D.M. Spink % Copyright (c) 2000. if iscell(crv1.knots) | iscell(crv2.knots) error('Both NURBS must be curves'); end % ensure both curves have a common degree d = max([crv1.order, crv2.order]); crv1 = nrbdegelev(crv1, d - crv1.order); crv2 = nrbdegelev(crv2, d - crv2.order); % merge the knot vectors, to obtain a common knot vector k1 = crv1.knots; k2 = crv2.knots; ku = unique([k1 k2]); n = length(ku); ka = []; kb = []; for i = 1:n i1 = length(find(k1 == ku(i))); i2 = length(find(k2 == ku(i))); m = max(i1, i2); ka = [ka ku(i)*ones(1,m-i1)]; kb = [kb ku(i)*ones(1,m-i2)]; end crv1 = nrbkntins(crv1, ka); crv2 = nrbkntins(crv2, kb); coefs(:,:,1) = crv1.coefs; coefs(:,:,2) = crv2.coefs; srf = nrbmak(coefs, {crv1.knots [0 0 1 1]});

nurbs_toolbox/nrbtestcrv.m

function crv = nrbtestcrv % Constructs a simple test curve. % D.M. Spink % Copyright (c) 2000 pnts = [0.5 1.5 4.5 3.0 7.5 6.0 8.5; 3.0 5.5 5.5 1.5 1.5 4.0 4.5; 0.0 0.0 0.0 0.0 0.0 0.0 0.0]; crv = nrbmak(pnts,[0 0 0 1/4 1/2 3/4 3/4 1 1 1]);

nurbs_toolbox/nrbtestsrf.m

function srf = nrbtestsrf % Constructs a simple test surface. % D.M. Spink % Copyright (c) 2000. % allocate multi-dimensional array of control points pnts = zeros(3,5,5); % define a grid of control points % in this case a regular grid of u,v points % pnts(3,u,v) % pnts(:,:,1) = [ 0.0 3.0 5.0 8.0 10.0; % w*x 0.0 0.0 0.0 0.0 0.0; % w*y 2.0 2.0 7.0 7.0 8.0]; % w*z pnts(:,:,2) = [ 0.0 3.0 5.0 8.0 10.0; 3.0 3.0 3.0 3.0 3.0; 0.0 0.0 5.0 5.0 7.0]; pnts(:,:,3) = [ 0.0 3.0 5.0 8.0 10.0; 5.0 5.0 5.0 5.0 5.0; 0.0 0.0 5.0 5.0 7.0]; pnts(:,:,4) = [ 0.0 3.0 5.0 8.0 10.0; 8.0 8.0 8.0 8.0 8.0; 5.0 5.0 8.0 8.0 10.0]; pnts(:,:,5) = [ 0.0 3.0 5.0 8.0 10.0; 10.0 10.0 10.0 10.0 10.0; 5.0 5.0 8.0 8.0 10.0]; % knots knots{1} = [0 0 0 1/3 2/3 1 1 1]; % knots along u knots{2} = [0 0 0 1/3 2/3 1 1 1]; % knots along v % make and draw nurbs surface srf = nrbmak(pnts,knots);

nurbs_toolbox/nrbtform.m

function nurbs = nrbtform(nurbs,tmat) % % Function Name: % % nrbtform - Apply transformation matrix to the NURBS. % % Calling Sequence: % % tnurbs = nrbtform(nurbs,tmatrix); % % Parameters: % % nurbs : NURBS data structure (see nrbmak for details). % % tmatrix : Transformation matrix, a matrix of size (4,4) defining % a single or multiple transformations. % % tnurbs : The return transformed NURBS data structure. % % Description: % % The NURBS is transform as defined a transformation matrix of size (4,4), % such as a rotation, translation or change in scale. The transformation % matrix can define a single transformation or multiple series of % transformations. The matrix can be simple constructed by the functions % vecscale, vectrans, vecrotx, vecroty, and vecrotz. % % Examples: % % Rotate a square by 45 degrees about the z axis. % % rsqr = nrbtform(nrbrect(), vecrotz(deg2rad(45))); % nrbplot(rsqr, 10); % % See: % % vecscale, vectrans, vecrotx, vecroty, vecrotz % D.M. Spink % Copyright (c) 2000 if nargin < 2 error('Not enough input arguments!'); end; if iscell(nurbs.knots) % NURBS is a surface [dim,nu,nv] = size(nurbs.coefs); nurbs.coefs = reshape(tmat*reshape(nurbs.coefs,dim,nu*nv),[dim nu nv]); else % NURBS is a curve nurbs.coefs = tmat*nurbs.coefs; end

nurbs_toolbox/nrbtransp.m

function tsrf = nrbtransp(srf) % % Function Name: % % nrbtransp - Transpose a NURBS surface, by swapping U and V directions. % % Calling Sequence: % % tsrf = nrbtransp(srf) % % Parameters: % % srf : NURBS surface, see nrbmak. % % tsrf : NURBS surface with U and V diretions transposed. % % Description: % % Utility function that transposes a NURBS surface, by swapping U and % V directions. NURBS curves cannot be transposed. % D.M. Spink % Copyright (c) 2000. if ~iscell(srf.knots) error(' A NURBSs curve cannot be transposed.'); end tsrf = nrbmak(permute(srf.coefs,[1 3 2]), fliplr(srf.knots));

nurbs_toolbox/rad2deg.m

function deg = rad2deg(rad) % % Function Name: % % rad2deg - Convert radians to degrees. % % Calling Sequence: % % rad = rad2deg(deg); % % Parameters: % % rad : Angle in radians. % % deg : Angle in degrees. % % Description: % % Convenient utility function for converting radians to degrees, which are % often the required angular units for functions in the NURBS toolbox. % % Examples: % % Convert 0.3 radians to degrees % % rad = deg2rad(0.3); % D.M. Spink % Copyright (c) 2000. deg = 180.0*rad/pi;

nurbs_toolbox/vecangle.m

function ang = vecangle(num,den) % % Function Name: % % vecangle - An alternative to atan, returning an arctangent in the % range 0 to 2*pi. % % Calling Sequence: % % ang = vecmag2(num,dum) % % Parameters: % % num : Numerator, vector of size (1,nv). % dem : Denominator, vector of size (1,nv). % ang : Arctangents, row vector of angles. % % Description: % % The components of the vector ang are the arctangent of the corresponding % enties of num./dem. This function is an alternative for % atan, returning an angle in the range 0 to 2*pi. % % Examples: % % Find the atan(1.2,2.0) and atan(1.5,3.4) using vecangle % % ang = vecangle([1.2 1.5], [2.0 3.4]); % D.M. Spink % Copyright (c) 2000. ang = atan2(num,den); index = find(ang < 0.0); ang(index) = 2*pi+ang(index);