Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[5628] trunk/octave-forge/extra/nurbs
|
From: <cd...@us...> - 2009-03-22 13:20:58
|
Revision: 5628
http://octave.svn.sourceforge.net/octave/?rev=5628&view=rev
Author: cdf
Date: 2009-03-22 13:20:48 +0000 (Sun, 22 Mar 2009)
Log Message:
-----------
bspeval
Modified Paths:
--------------
trunk/octave-forge/extra/nurbs/INDEX
trunk/octave-forge/extra/nurbs/inst/bspdegelev.m
trunk/octave-forge/extra/nurbs/inst/findspan.m
Added Paths:
-----------
trunk/octave-forge/extra/nurbs/inst/nrbsrfgradient.m
trunk/octave-forge/extra/nurbs/src/
trunk/octave-forge/extra/nurbs/src/Makefile
trunk/octave-forge/extra/nurbs/src/basisfun.c
trunk/octave-forge/extra/nurbs/src/bspdegelev.c
trunk/octave-forge/extra/nurbs/src/bspderiv.c
trunk/octave-forge/extra/nurbs/src/bspeval.c
trunk/octave-forge/extra/nurbs/src/bspkntins.c
trunk/octave-forge/extra/nurbs/src/findspan.c
trunk/octave-forge/extra/nurbs/src/mexmat.c
trunk/octave-forge/extra/nurbs/src/mexmat.h
Modified: trunk/octave-forge/extra/nurbs/INDEX
===================================================================
--- trunk/octave-forge/extra/nurbs/INDEX 2009-03-22 12:36:28 UTC (rev 5627)
+++ trunk/octave-forge/extra/nurbs/INDEX 2009-03-22 13:20:48 UTC (rev 5628)
@@ -19,6 +19,7 @@
nrbruled
nrbcoons
nrbplot
+ nrbsrfgradient
Low level functions
bspeval
bspderiv
Modified: trunk/octave-forge/extra/nurbs/inst/bspdegelev.m
===================================================================
--- trunk/octave-forge/extra/nurbs/inst/bspdegelev.m 2009-03-22 12:36:28 UTC (rev 5627)
+++ trunk/octave-forge/extra/nurbs/inst/bspdegelev.m 2009-03-22 13:20:48 UTC (rev 5628)
@@ -13,53 +13,284 @@
## 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
+function [ic,ik] = bspdegelev(d,c,k,t)
+% BSPDEGELEV Degree elevate a univariate B-Spline.
% -------------------------------------------------------------------------
-% ADAPTATION of BASISFUN from C Routine
+% ADAPTATION of BSPDEGELEV 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;
+% [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.
+%
- 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 % }
+[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; %
- N(j+1) = saved; % N[j] = saved;
-end % }
-
- % mxFree(left);
- % mxFree(right);
- % }
+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
Modified: trunk/octave-forge/extra/nurbs/inst/findspan.m
===================================================================
--- trunk/octave-forge/extra/nurbs/inst/findspan.m 2009-03-22 12:36:28 UTC (rev 5627)
+++ trunk/octave-forge/extra/nurbs/inst/findspan.m 2009-03-22 13:20:48 UTC (rev 5628)
@@ -1,44 +1,47 @@
-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);
- % }
+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
+%
+% U(1) <= u <= U(end)
+% RETURN:
+%
+% s - knot span
+%
+% Algorithm A2.1 from 'The NURBS BOOK' pg68
+
+ % int findspan(int n, int p, double u, double *U) {
+if ((nargin ~= 4) || (u<U(1)) || (u>U(end)))
+ print_usage ()
+end
+ % 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);
+ % }
Added: trunk/octave-forge/extra/nurbs/inst/nrbsrfgradient.m
===================================================================
--- trunk/octave-forge/extra/nurbs/inst/nrbsrfgradient.m (rev 0)
+++ trunk/octave-forge/extra/nurbs/inst/nrbsrfgradient.m 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,49 @@
+## Copyright (C) 2009 Carlo de Falco
+##
+## 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, write to the Free Software
+## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{dzdx}, @var{dzdy}, @var{z}, @var{x}, @var{y}]=} nrbsrfgradient (@var{nrb}, @var{nrbder}, @var{u}, @var{v})
+## Compute the gradient of a NURBS surface.
+## @seealso{nrbderiv}
+## @end deftypefn
+
+## Author: Carlo de Falco <carlo@guglielmo.local>
+## Created: 2009-03-17
+
+function [dzdx, dzdy] = nrbsrfgradient (nrb, nrbder, u, v)
+
+ if ((nargin != 4) || (nargout>2))
+ print_usage();
+ endif
+
+ [np, dp] = nrbdeval (nrb, nrbder, {u, v});
+
+ dxdu = squeeze (dp{1}(1,:,:));
+ dydu = squeeze (dp{1}(2,:,:));
+ dzdu = squeeze (dp{1}(3,:,:));
+ dxdv = squeeze (dp{2}(1,:,:));
+ dydv = squeeze (dp{2}(2,:,:));
+ dzdv = squeeze (dp{2}(3,:,:));
+
+ dzdx = dzdy = zeros(length(u), length(v));
+
+ dzdx(dxdu~=0) += (dzdu./dxdu)(dxdu~=0);
+ dzdx(dxdv~=0) += (dzdv./dxdv)(dxdv~=0);
+
+ dzdy(dydu~=0) += (dzdu./dydu)(dydu~=0);
+ dzdy(dydv~=0) += (dzdv./dydv)(dydv~=0);
+
+endfunction
Added: trunk/octave-forge/extra/nurbs/src/Makefile
===================================================================
--- trunk/octave-forge/extra/nurbs/src/Makefile (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/Makefile 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,20 @@
+mex=mkoctfile --mex
+
+all: bspeval.mex bspderiv.mex bspkntins.mex bspdegelev.mex
+
+bspeval.mex: bspeval.c basisfun.c findspan.c mexmat.c
+ $(mex) bspeval.c basisfun.c findspan.c mexmat.c
+
+bspderiv.mex: bspderiv.c mexmat.c
+ $(mex) bspderiv.c mexmat.c
+
+bspkntins.mex: bspkntins.c findspan.c mexmat.c
+ $(mex) bspkntins.c findspan.c mexmat.c
+
+bspdegelev.mex: bspdegelev.c mexmat.c
+ $(mex) bspdegelev.c mexmat.c
+
+clean:
+ rm -f bspeval.mex bspeval.o bspderiv.mex \
+ bspderiv.o bspkntins.mex bspkntins.o bspdegelev.mex bspdegelev.o
+
Property changes on: trunk/octave-forge/extra/nurbs/src/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/basisfun.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/basisfun.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/basisfun.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,46 @@
+// Basis Function.
+//
+// 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.
+
+#include "mexmat.h"
+
+void basisfun(int i, double u, int p, double *U, double *N)
+{
+ int j,r;
+ double saved, temp;
+
+ // work space
+ double *left = (double*) mxMalloc((p+1)*sizeof(double));
+ double *right = (double*) mxMalloc((p+1)*sizeof(double));
+
+ N[0] = 1.0;
+ for (j = 1; j <= p; j++)
+ {
+ left[j] = u - U[i+1-j];
+ right[j] = U[i+j] - u;
+ saved = 0.0;
+
+ for (r = 0; r < j; r++)
+ {
+ temp = N[r] / (right[r+1] + left[j-r]);
+ N[r] = saved + right[r+1] * temp;
+ saved = left[j-r] * temp;
+ }
+
+ N[j] = saved;
+ }
+
+ mxFree(left);
+ mxFree(right);
+}
Property changes on: trunk/octave-forge/extra/nurbs/src/basisfun.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/bspdegelev.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/bspdegelev.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/bspdegelev.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,349 @@
+
+#include "mexmat.h"
+
+#define min(x,y) (x < y ? x : y)
+#define max(x,y) (x > y ? x : y)
+
+// Compute logarithm of the gamma function
+// Algorithm from 'Numerical Recipes in C, 2nd Edition' pg214.
+double gammaln(double xx)
+{
+ double x,y,tmp,ser;
+ static double cof[6] = {76.18009172947146,-86.50532032291677,
+ 24.01409824083091,-1.231739572450155,
+ 0.12086650973866179e-2, -0.5395239384953e-5};
+ int j;
+ y = x = xx;
+ tmp = x + 5.5;
+ tmp -= (x+0.5) * log(tmp);
+ ser = 1.000000000190015;
+ for (j=0; j<=5; j++) ser += cof[j]/++y;
+ return -tmp+log(2.5066282746310005*ser/x);
+}
+
+// computes ln(n!)
+// Numerical Recipes in C
+// Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215.
+double factln(int n)
+{
+ static int ntop = 0;
+ static double a[101];
+
+ if (n <= 1) return 0.0;
+ while (n > ntop)
+ {
+ ++ntop;
+ a[ntop] = gammaln(ntop+1.0);
+ }
+ return a[n];
+}
+
+// Computes the binomial coefficient.
+//
+// ( n ) n!
+// ( ) = --------
+// ( k ) k!(n-k)!
+//
+// Algorithm from 'Numerical Recipes in C, 2nd Edition' pg215.
+double bincoeff(int n, int k)
+{
+ return floor(0.5+exp(factln(n)-factln(k)-factln(n-k)));
+}
+
+// Degree elevate a B-Spline t times.
+// i.e. from d to d+t
+//
+// INPUT:
+//
+// n,p,U,Pw,t
+//
+// OUTPUT:
+//
+// nh,Uh,Qw
+//
+// Modified version of Algorithm A5.9 from 'The NURBS BOOK' pg206.
+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);
+ double **ictrl = vec2mat(ic, mc, nc*(t+1));
+
+ n = nc - 1;
+
+ bezalfs = matrix(d+1,d+t+1);
+ bpts = matrix(mc,d+1);
+ ebpts = matrix(mc,d+t+1);
+ Nextbpts = matrix(mc,d+1);
+ alfs = (double *) mxMalloc(d*sizeof(double));
+
+ m = n + d + 1;
+ ph = d + t;
+ ph2 = ph / 2;
+
+ // compute bezier degree elevation coefficeients
+ bezalfs[0][0] = bezalfs[ph][d] = 1.0;
+
+ for (i = 1; i <= ph2; i++)
+ {
+ inv = 1.0 / bincoeff(ph,i);
+ mpi = min(d,i);
+
+ for (j = max(0,i-t); j <= mpi; j++)
+ bezalfs[i][j] = inv * bincoeff(d,j) * bincoeff(t,i-j);
+ }
+
+ for (i = ph2+1; i <= ph-1; i++)
+ {
+ mpi = min(d, i);
+ for (j = max(0,i-t); j <= mpi; j++)
+ bezalfs[i][j] = bezalfs[ph-i][d-j];
+ }
+
+ mh = ph;
+ kind = ph+1;
+ r = -1;
+ a = d;
+ b = d+1;
+ cind = 1;
+ ua = k[0];
+
+ for (ii = 0; ii < mc; ii++)
+ ictrl[0][ii] = ctrl[0][ii];
+
+ for (i = 0; i <= ph; i++)
+ ik[i] = ua;
+
+ // initialise first bezier seg
+ for (i = 0; i <= d; i++)
+ for (ii = 0; ii < mc; ii++)
+ bpts[i][ii] = ctrl[i][ii];
+
+ // big loop thru knot vector
+ while (b < m)
+ {
+ i = b;
+ while (b < m && k[b] == k[b+1])
+ b++;
+
+ mul = b - i + 1;
+ mh += mul + t;
+ ub = k[b];
+ oldr = r;
+ r = d - mul;
+
+ // insert knot u(b) r times
+ if (oldr > 0)
+ lbz = (oldr+2) / 2;
+ else
+ lbz = 1;
+
+ if (r > 0)
+ rbz = ph - (r+1)/2;
+ else
+ rbz = ph;
+
+ if (r > 0)
+ {
+ // insert knot to get bezier segment
+ numer = ub - ua;
+ for (q = d; q > mul; q--)
+ alfs[q-mul-1] = numer / (k[a+q]-ua);
+ for (j = 1; j <= r; j++)
+ {
+ save = r - j;
+ s = mul + j;
+
+ for (q = d; q >= s; q--)
+ for (ii = 0; ii < mc; ii++)
+ bpts[q][ii] = alfs[q-s]*bpts[q][ii]+(1.0-alfs[q-s])*bpts[q-1][ii];
+
+ for (ii = 0; ii < mc; ii++)
+ Nextbpts[save][ii] = bpts[d][ii];
+ }
+ }
+ // end of insert knot
+
+ // degree elevate bezier
+ for (i = lbz; i <= ph; i++)
+ {
+ for (ii = 0; ii < mc; ii++)
+ ebpts[i][ii] = 0.0;
+ mpi = min(d, i);
+ for (j = max(0,i-t); j <= mpi; j++)
+ for (ii = 0; ii < mc; ii++)
+ ebpts[i][ii] = ebpts[i][ii] + bezalfs[i][j]*bpts[j][ii];
+ }
+ // end of degree elevating bezier
+
+ if (oldr > 1)
+ {
+ // must remove knot u=k[a] oldr times
+ first = kind - 2;
+ last = kind;
+ den = ub - ua;
+ bet = (ub-ik[kind-1]) / den;
+
+ // knot removal loop
+ for (tr = 1; tr < oldr; tr++)
+ {
+ i = first;
+ j = last;
+ kj = j - kind + 1;
+ while (j - i > tr)
+ {
+ // loop and compute the new control points
+ // for one removal step
+ if (i < cind)
+ {
+ alf = (ub-ik[i])/(ua-ik[i]);
+ for (ii = 0; ii < mc; ii++)
+ ictrl[i][ii] = alf * ictrl[i][ii] + (1.0-alf) * ictrl[i-1][ii];
+ }
+ if (j >= lbz)
+ {
+ if (j-tr <= kind-ph+oldr)
+ {
+ gam = (ub-ik[j-tr]) / den;
+ for (ii = 0; ii < mc; ii++)
+ ebpts[kj][ii] = gam*ebpts[kj][ii] + (1.0-gam)*ebpts[kj+1][ii];
+ }
+ else
+ {
+ for (ii = 0; ii < mc; ii++)
+ ebpts[kj][ii] = bet*ebpts[kj][ii] + (1.0-bet)*ebpts[kj+1][ii];
+ }
+ }
+ i++;
+ j--;
+ kj--;
+ }
+
+ first--;
+ last++;
+ }
+ }
+ // end of removing knot n=k[a]
+
+ // load the knot ua
+ if (a != d)
+ for (i = 0; i < ph-oldr; i++)
+ {
+ ik[kind] = ua;
+ kind++;
+ }
+
+ // load ctrl pts into ic
+ for (j = lbz; j <= rbz; j++)
+ {
+ for (ii = 0; ii < mc; ii++)
+ ictrl[cind][ii] = ebpts[j][ii];
+ cind++;
+ }
+
+ if (b < m)
+ {
+ // setup for next pass thru loop
+ for (j = 0; j < r; j++)
+ for (ii = 0; ii < mc; ii++)
+ bpts[j][ii] = Nextbpts[j][ii];
+ for (j = r; j <= d; j++)
+ for (ii = 0; ii < mc; ii++)
+ bpts[j][ii] = ctrl[b-d+j][ii];
+ a = b;
+ b++;
+ ua = ub;
+ }
+ else
+ // end knot
+ for (i = 0; i <= ph; i++)
+ ik[kind+i] = ub;
+ }
+ // end while loop
+
+ *nh = mh - ph - 1;
+
+ freevec2mat(ctrl);
+ freevec2mat(ictrl);
+ freematrix(bezalfs);
+ freematrix(bpts);
+ freematrix(ebpts);
+ freematrix(Nextbpts);
+ mxFree(alfs);
+
+ return(ierr);
+}
+
+// Matlab gateway function
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+ int mc, nc;
+ int mk, nk;
+ int mu, nu;
+
+ int d, t;
+ double *c, *k;
+
+ int nic, nik;
+ double *ic, *ik;
+
+ int nh, i, n;
+ double *inc, *ink;
+
+ if (nrhs != 4)
+ mexErrMsgTxt("Four inputs required\n");
+
+ if (nlhs > 2)
+ mexErrMsgTxt("Maximum of two outputs required\n");
+
+ // spline degree
+ d = (int) mxGetScalar(prhs[0]);
+
+ // control points
+ c = mxGetPr(prhs[1]);
+ mc = mxGetM(prhs[1]);
+ nc = mxGetN(prhs[1]);
+
+ // knot sequence
+ k = mxGetPr(prhs[2]);
+ mk = mxGetM(prhs[2]);
+ nk = mxGetN(prhs[2]);
+
+ // raise degree t times
+ t = (int) mxGetScalar(prhs[3]);
+
+ // allocate work space t times larger than original number
+ // of control points and knots
+ inc = (double *) mxMalloc((t+1)*mc*nc*sizeof(double));
+ ink = (double *) mxMalloc((t+1)*nk*sizeof(double));
+
+ bspdegelev(d, c, mc, nc, k, nk, t, &nh, inc, ink);
+
+ // new coefficients and knot sequence
+ nic = nh + 1;
+ nik = nic + d + t + 1;
+ plhs[0] = mxCreateDoubleMatrix(mc, nic, mxREAL);
+ ic = mxGetPr(plhs[0]);
+ plhs[1] = mxCreateDoubleMatrix(mk, nik, mxREAL);
+ ik = mxGetPr(plhs[1]);
+
+ // copy new control points
+ n = mc * nic;
+ for (i = 0; i < n; i++)
+ ic[i] = inc[i];
+
+ // copy new knots
+ for (i = 0; i < nik; i++)
+ ik[i] = ink[i];
+
+ // free working space
+ mxFree(ink);
+ mxFree(inc);
+}
Property changes on: trunk/octave-forge/extra/nurbs/src/bspdegelev.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/bspderiv.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/bspderiv.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/bspderiv.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,94 @@
+// Compute the control points and knot sequence a univariate B-Spline
+// derivative.
+//
+// MATLAB SYNTAX:
+//
+// [dc,dk] = bspderiv(d,c,k)
+//
+// INPUT:
+//
+// d - degree of the B-Spline
+// c - control points double matrix(mc,nc)
+// k - knot sequence double vector(nk)
+//
+// OUTPUT:
+//
+// dc - control points of the derivative double matrix(mc,nc)
+// dk - knot sequence of the derivative double vector(nk)
+//
+//
+
+#include "mexmat.h"
+
+// Modified version of Algorithm A3.3 from 'The NURBS BOOK' pg98.
+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
+ double **dctrl = vec2mat(dc,mc,nc-1);
+
+ for (i = 0; i < nc-1; i++) {
+ tmp = d / (k[i+d+1] - k[i+1]);
+ for (j = 0; j < mc; j++) {
+ dctrl[i][j] = tmp * (ctrl[i+1][j] - ctrl[i][j]);
+ }
+ }
+
+ j = 0;
+ for (i = 1; i < nk-1; i++)
+ dk[j++] = k[i];
+
+ freevec2mat(dctrl);
+ freevec2mat(ctrl);
+
+ return ierr;
+}
+
+// Matlab gateway function
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+ int mc, nc;
+ int mk, nk;
+ int mu, nu;
+
+ int d;
+ double *c, *p, *k, *dc, *dk;
+
+ if (nrhs != 3)
+ mexErrMsgTxt("Three inputs required\n");
+
+ if (nlhs > 2)
+ mexErrMsgTxt("Two output required maximum\n");
+
+ // spline degree
+ d = (int) mxGetScalar(prhs[0]);
+
+ // control points
+ c = mxGetPr(prhs[1]);
+ mc = mxGetM(prhs[1]);
+ nc = mxGetN(prhs[1]);
+
+ // knot sequence
+ k = mxGetPr(prhs[2]);
+ mk = mxGetM(prhs[2]);
+ nk = mxGetN(prhs[2]);
+
+ // control points of the derivative
+ plhs[0] = mxCreateDoubleMatrix(mc, nc-1, mxREAL);
+ dc = mxGetPr(plhs[0]);
+
+ // knot sequence of the derivative
+ plhs[1] = mxCreateDoubleMatrix(mk, nk-2, mxREAL);
+ dk = mxGetPr(plhs[1]);
+
+ bspderiv(d, c, mc, nc, k, nk, dc, dk);
+}
+
+
+
Property changes on: trunk/octave-forge/extra/nurbs/src/bspderiv.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/bspeval.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/bspeval.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/bspeval.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,105 @@
+// Evaluation of univariate B-Spline.
+//
+// MATLAB SYNTAX:
+//
+// p = bspeval(d,c,k,u)
+//
+// INPUT:
+//
+// d - degree of the B-Spline integer
+// c - control points double matrix(mc,nc)
+// k - knot sequence double vector(nk)
+// u - parametric points double vector(nu)
+//
+// OUTPUT:
+//
+// p - evaluated points double matrix(mc,nu)
+//
+//
+
+#include "mexmat.h"
+
+// Modified version of Algorithm A3.1 from 'The NURBS BOOK' pg82.
+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
+ double **pnt = vec2mat(p,mc,nu);
+
+ // space for the basis functions
+ double *N = (double*) mxMalloc((d+1)*sizeof(double));
+
+ // for each parametric point i
+ for (col = 0; col < nu; col++)
+ {
+ // find the span of u[col]
+ s = findspan(nc-1, d, u[col], k);
+ basisfun(s, u[col], d, k, N);
+
+ tmp1 = s - d;
+ for (row = 0; row < mc; row++)
+ {
+ tmp2 = 0.0;
+ for (i = 0; i <= d; i++)
+ tmp2 += N[i] * ctrl[tmp1+i][row];
+ pnt[col][row] = tmp2;
+ }
+ }
+
+ mxFree(N);
+ freevec2mat(pnt);
+ freevec2mat(ctrl);
+
+ return ierr;
+}
+
+// Matlab gateway function
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+ int mc, nc;
+ int mk, nk;
+ int mu, nu;
+
+ int d;
+ double *c, *p, *k, *u;
+
+ if (nrhs != 4)
+ mexErrMsgTxt("Four inputs required\n");
+
+ if (nlhs > 1)
+ mexErrMsgTxt("One output required\n");
+
+ // spline degree
+ d = (int) mxGetScalar(prhs[0]);
+
+ // control points
+ c = mxGetPr(prhs[1]);
+ mc = mxGetM(prhs[1]);
+ nc = mxGetN(prhs[1]);
+
+ // knot sequence
+ k = mxGetPr(prhs[2]);
+ mk = mxGetM(prhs[2]);
+ nk = mxGetN(prhs[2]);
+
+ // parametric points
+ u = mxGetPr(prhs[3]);
+ mu = mxGetM(prhs[3]);
+ nu = mxGetN(prhs[3]);
+
+ // evaluated points
+ plhs[0] = mxCreateDoubleMatrix(mc, nu, mxREAL);
+ p = mxGetPr(plhs[0]);
+
+ bspeval(d, c, mc, nc, k, nk, u, nu, p);
+}
+
+
+
Property changes on: trunk/octave-forge/extra/nurbs/src/bspeval.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/bspkntins.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/bspkntins.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/bspkntins.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,130 @@
+// Insert Knot into a B-Spline
+//
+// INPUT:
+//
+// d - spline degree integer
+// c - control points double matrix(mc,nc)
+// k - knot sequence double vector(nk)
+// u - new knots double vector(nu)
+//
+// OUTPUT:
+//
+// ic - new control points double matrix(mc,nc+nu)
+// ik - new knot sequence double vector(nk+nu)
+//
+// Modified version of Algorithm A5.4 from 'The NURBS BOOK' pg164.
+
+#include "mexmat.h"
+
+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);
+ double **ictrl = vec2mat(ic, mc, nc+nu);
+
+ n = nc - 1;
+ r = nu - 1;
+
+ m = n + d + 1;
+ a = findspan(n, d, u[0], k);
+ b = findspan(n, d, u[r], k);
+ ++b;
+
+ for (q = 0; q < mc; q++)
+ {
+ for (j = 0; j <= a-d; j++) ictrl[j][q] = ctrl[j][q];
+ for (j = b-1; j <= n; j++) ictrl[j+r+1][q] = ctrl[j][q];
+ }
+ for (j = 0; j <= a; j++) ik[j] = k[j];
+ for (j = b+d; j <= m; j++) ik[j+r+1] = k[j];
+
+ i = b + d - 1;
+ s = b + d + r;
+ for (j = r; j >= 0; j--)
+ {
+ while (u[j] <= k[i] && i > a)
+ {
+ for (q = 0; q < mc; q++)
+ ictrl[s-d-1][q] = ctrl[i-d-1][q];
+ ik[s] = k[i];
+ --s;
+ --i;
+ }
+ for (q = 0; q < mc; q++)
+ ictrl[s-d-1][q] = ictrl[s-d][q];
+ for (l = 1; l <= d; l++)
+ {
+ ind = s - d + l;
+ alfa = ik[s+l] - u[j];
+ if (fabs(alfa) == 0.0)
+ for (q = 0; q < mc; q++)
+ ictrl[ind-1][q] = ictrl[ind][q];
+ else
+ {
+ alfa /= (ik[s+l] - k[i-d+l]);
+ for (q = 0; q < mc; q++)
+ ictrl[ind-1][q] = alfa*ictrl[ind-1][q]+(1.0-alfa)*ictrl[ind][q];
+ }
+ }
+
+ ik[s] = u[j];
+ --s;
+ }
+
+ freevec2mat(ctrl);
+ freevec2mat(ictrl);
+
+ return ierr;
+}
+
+// Matlab gateway function
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+ int mc, nc;
+ int mk, nk;
+ int mu, nu;
+
+ int d;
+ double *c, *k, *u;
+
+ int nic, nik;
+ double *ic, *ik;
+
+ if (nrhs != 4)
+ mexErrMsgTxt("Four inputs required\n");
+
+ if (nlhs > 2)
+ mexErrMsgTxt("Maximum of two outputs required\n");
+
+ // spline degree
+ d = (int) mxGetScalar(prhs[0]);
+
+ // control points
+ c = mxGetPr(prhs[1]);
+ mc = mxGetM(prhs[1]);
+ nc = mxGetN(prhs[1]);
+
+ // knot sequence
+ k = mxGetPr(prhs[2]);
+ mk = mxGetM(prhs[2]);
+ nk = mxGetN(prhs[2]);
+
+ // knots to be inserted
+ u = mxGetPr(prhs[3]);
+ mu = mxGetM(prhs[3]);
+ nu = mxGetN(prhs[3]);
+
+ // new coefficients and knot sequence
+ nic = nc + nu;
+ nik = nk + nu;
+ plhs[0] = mxCreateDoubleMatrix(mc, nic, mxREAL);
+ ic = mxGetPr(plhs[0]);
+ plhs[1] = mxCreateDoubleMatrix(mk, nik, mxREAL);
+ ik = mxGetPr(plhs[1]);
+
+ bspkntins(d, c, mc, nc, k, nk, u, nu, ic, ik);
+}
Property changes on: trunk/octave-forge/extra/nurbs/src/bspkntins.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/findspan.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/findspan.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/findspan.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,38 @@
+// Find the knot span of the parametric point 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+1]) return(n);
+
+ // do binary search
+ low = p;
+ high = n + 1;
+ mid = (low + high) / 2;
+ while (u < U[mid] || u >= U[mid+1])
+ {
+ if (u < U[mid])
+ high = mid;
+ else
+ low = mid;
+ mid = (low + high) / 2;
+ }
+
+ return(mid);
+}
+
Property changes on: trunk/octave-forge/extra/nurbs/src/findspan.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/mexmat.c
===================================================================
--- trunk/octave-forge/extra/nurbs/src/mexmat.c (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/mexmat.c 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,43 @@
+// Useful functions for accessing and manipulating matlab data structures.
+// =======================================================================
+
+#include "mexmat.h"
+
+// convert c vector to c matrix.
+double **vec2mat(double *vec, int nrows, int ncols)
+{
+ int col;
+ double **mat;
+
+ mat = (double**) mxMalloc (ncols*sizeof(double*));
+ mat[0] = vec;
+ for (col = 1; col < ncols; col++)
+ mat[col] = mat[col-1] + nrows;
+ return mat;
+}
+
+// create a new c matrix
+double **matrix(int nrows, int ncols)
+{
+ int col;
+ double **mat;
+
+ mat = (double**) mxMalloc (ncols*sizeof(double*));
+ mat[0] = (double*) mxMalloc (nrows*ncols*sizeof(double));
+ for (col = 1; col < ncols; col++)
+ mat[col] = mat[col-1] + nrows;
+ return mat;
+}
+
+void freevec2mat(double **mat)
+{
+ mxFree(mat);
+}
+
+void freematrix(double **mat)
+{
+ mxFree(mat[0]);
+ mxFree(mat);
+}
+
+
Property changes on: trunk/octave-forge/extra/nurbs/src/mexmat.c
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/extra/nurbs/src/mexmat.h
===================================================================
--- trunk/octave-forge/extra/nurbs/src/mexmat.h (rev 0)
+++ trunk/octave-forge/extra/nurbs/src/mexmat.h 2009-03-22 13:20:48 UTC (rev 5628)
@@ -0,0 +1,16 @@
+// Useful functions for accessing and manipulating matlab data structures.
+// =======================================================================
+
+#ifndef __MEXMAT_H
+#define __MEXMAT_H
+
+#include <math.h>
+#include "mex.h"
+
+double **vec2mat(double *vec, int nrows, int ncols);
+double **matrix(int nrows, int ncols);
+void freevec2mat(double **mat);
+void freematrix(double **mat);
+
+#endif // __MEXMAT_H
+
Property changes on: trunk/octave-forge/extra/nurbs/src/mexmat.h
___________________________________________________________________
Added: svn:executable
+ *
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|