From: <cd...@us...> - 2009-03-16 14:01:14
|
Revision: 5612 http://octave.svn.sourceforge.net/octave/?rev=5612&view=rev Author: cdf Date: 2009-03-16 13:37:04 +0000 (Mon, 16 Mar 2009) Log Message: ----------- removed calls to isstr Modified Paths: -------------- trunk/octave-forge/extra/nurbs/inst/bspdegelev.m trunk/octave-forge/extra/nurbs/inst/bspderiv.m trunk/octave-forge/extra/nurbs/inst/nrbdemo.m trunk/octave-forge/extra/nurbs/inst/nrbplot.m Modified: trunk/octave-forge/extra/nurbs/inst/bspdegelev.m =================================================================== --- trunk/octave-forge/extra/nurbs/inst/bspdegelev.m 2009-03-16 07:16:40 UTC (rev 5611) +++ trunk/octave-forge/extra/nurbs/inst/bspdegelev.m 2009-03-16 13:37:04 UTC (rev 5612) @@ -1,296 +1,65 @@ -## Copyright (C) 2003 Mark Spink, 2007 Daniel Claxton -## -## This program is free software; you can redistribute it and/or modify -## it under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 2 of the License, or -## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, -## but WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with this program; if not, see <http://www.gnu.org/licenses/>. - -function [ic,ik] = bspdegelev(d,c,k,t) -% BSPDEGELEV Degree elevate a univariate B-Spline. -% ------------------------------------------------------------------------- -% ADAPTATION of BSPDEGELEV from C Routine -% ------------------------------------------------------------------------- -% -% 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. -% - -[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)); \ No newline at end of file +## Copyright (C) 2003 Mark Spink, 2007 Daniel Claxton +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, see <http://www.gnu.org/licenses/>. + +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); + % } Modified: trunk/octave-forge/extra/nurbs/inst/bspderiv.m =================================================================== --- trunk/octave-forge/extra/nurbs/inst/bspderiv.m 2009-03-16 07:16:40 UTC (rev 5611) +++ trunk/octave-forge/extra/nurbs/inst/bspderiv.m 2009-03-16 13:37:04 UTC (rev 5612) @@ -67,4 +67,4 @@ % freevec2mat(ctrl); % % return ierr; - % } \ No newline at end of file + % } Modified: trunk/octave-forge/extra/nurbs/inst/nrbdemo.m =================================================================== --- trunk/octave-forge/extra/nurbs/inst/nrbdemo.m 2009-03-16 07:16:40 UTC (rev 5611) +++ trunk/octave-forge/extra/nurbs/inst/nrbdemo.m 2009-03-16 13:37:04 UTC (rev 5612) @@ -31,10 +31,10 @@ 'Test Extrusion','demoextrude', ... 'Test Revolution - Cup','demorevolve', ... 'Test Revolution - Ball & Torus','demotorus', ... - '2D Geomtery','demogeom'} + '2D Geomtery','demogeom'}; for ii=1:2:length(lst) printf("Demo number %d: %s\n", (ii+1)/2, lst{ii}) eval(lst{ii+1}) pause -endfor \ No newline at end of file +endfor Modified: trunk/octave-forge/extra/nurbs/inst/nrbplot.m =================================================================== --- trunk/octave-forge/extra/nurbs/inst/nrbplot.m 2009-03-16 07:16:40 UTC (rev 5611) +++ trunk/octave-forge/extra/nurbs/inst/nrbplot.m 2009-03-16 13:37:04 UTC (rev 5612) @@ -68,7 +68,7 @@ for i=3:2:nargs Param = eval(['p' int2str((i-3)/2 +1)]); Value = eval(['v' int2str((i-3)/2 +1)]); - if ~isstr(Param) + if ~ischar(Param) error('Parameter must be a string') elseif size(Param,1)~=1 error('Parameter must be a non-empty single row string.') @@ -76,13 +76,13 @@ switch lower(Param) case 'light' light = lower(Value); - if ~isstr(light) + if ~ischar(light) error('light must be a string.') elseif ~(strcmp(light,'off') | strcmp(light,'on')) error('light must be off | on') end case 'colormap' - if isstr(Value) + if ischar(Value) cmap = lower(Value); elseif size(Value,2) ~= 3 error('colormap must be a string or have exactly three columns.') This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |