--- a/extra/nurbs/inst/nrbderiv.m	Fri Jul 22 13:22:03 2011 +0000
+++ b/extra/nurbs/inst/nrbderiv.m	Fri Jul 22 15:17:10 2011 +0000
@@ -1,11 +1,12 @@
-function dnurbs = nrbderiv(nurbs)
+function varargout = nrbderiv (nurbs)
 % 
-% NRBDERIV: Construct the first derivative representation of a
+% NRBDERIV2: Construct the first and second derivative representation of a
 %           NURBS curve, surface or volume.
 % 
 % Calling Sequence:
 % 
-%   ders = nrbderiv(nrb);
+%   ders = nrbderiv (nrb);
+%   [ders, ders2] = nrbderiv (nrb);
 % 
 % INPUT:
 % 
@@ -13,34 +14,28 @@
 %
 % OUTPUT:
 % 
-%   ders	: A data structure that represents the first
-% 		  derivatives of a NURBS curve, surface or volume.
+%   ders:  A data structure that represents the first
+% 		    derivatives of a NURBS curve, surface or volume.
+%   ders2: A data structure that represents the second
+% 		    derivatives of a NURBS curve, surface or volume.
 % 
 % 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.
+%   NURBS is more complex. We have followed in this function the same idea
+%   that was already used for the first derivative in the function nrbderiv.
+%   The second derivative data structure can be evaluated later with the
+%   function nrbdeval2.
 % 
 % See also:
 % 
 %       nrbdeval
 %
 % Copyright (C) 2000 Mark Spink
-% Copyright (C) 2010 Rafael Vazquez
 % Copyright (C) 2010 Carlo de Falco
+% Copyright (C) 2010, 2011 Rafael Vazquez
 %
 %    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
@@ -55,18 +50,24 @@
 %    You should have received a copy of the GNU General Public License
 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-if ~isstruct(nurbs)
+if (~isstruct(nurbs))
   error('NURBS representation is not structure!');
 end
 
-if ~strcmp(nurbs.form,'B-NURBS')
+if (~strcmp(nurbs.form,'B-NURBS'))
   error('Not a recognised NURBS representation');
 end
 
+% We raise the degree to avoid errors in the computation of the second derivative
+if (nargout == 2)
+  degelev  = max ([2 2] - (nurbs.order-1), 0);
+  nurbs    = nrbdegelev (nurbs, degelev);
+end
+
 degree = nurbs.order - 1;
 
-if iscell(nurbs.knots)
-  if size(nurbs.knots,2) == 3
+if (iscell(nurbs.knots))
+  if (size(nurbs.knots,2) == 3)
   % NURBS structure represents a volume
     num1 = nurbs.number(1);
     num2 = nurbs.number(2);
@@ -74,56 +75,98 @@
 
   % taking derivatives along the u direction
     dknots = nurbs.knots;
-    dcoefs = permute(nurbs.coefs,[1 3 4 2]);
-    dcoefs = reshape(dcoefs,4*num2*num3,num1);
-    [dcoefs,dknots{1}] = bspderiv(degree(1),dcoefs,nurbs.knots{1});
-    dcoefs = permute(reshape(dcoefs,[4 num2 num3 size(dcoefs,2)]),[1 4 2 3]);
-    dnurbs{1} = nrbmak(dcoefs, dknots);
+    dcoefs = permute (nurbs.coefs,[1 3 4 2]);
+    dcoefs = reshape (dcoefs,4*num2*num3,num1);
+    [dcoefs,dknots{1}] = bspderiv (degree(1),dcoefs,nurbs.knots{1});
+    dcoefs = permute (reshape (dcoefs,[4 num2 num3 size(dcoefs,2)]),[1 4 2 3]);
+    dnurbs{1} = nrbmak (dcoefs, dknots);
 
   % taking derivatives along the v direction
     dknots = nurbs.knots;
-    dcoefs = permute(nurbs.coefs,[1 2 4 3]);
-    dcoefs = reshape(dcoefs,4*num1*num3,num2);
-    [dcoefs,dknots{2}] = bspderiv(degree(2),dcoefs,nurbs.knots{2});
-    dcoefs = permute(reshape(dcoefs,[4 num1 num3 size(dcoefs,2)]),[1 2 4 3]);
-    dnurbs{2} = nrbmak(dcoefs, dknots);
+    dcoefs = permute (nurbs.coefs,[1 2 4 3]);
+    dcoefs = reshape (dcoefs,4*num1*num3,num2);
+    [dcoefs,dknots{2}] = bspderiv (degree(2),dcoefs,nurbs.knots{2});
+    dcoefs = permute (reshape (dcoefs,[4 num1 num3 size(dcoefs,2)]),[1 2 4 3]);
+    dnurbs{2} = nrbmak (dcoefs, dknots);
 
   % taking derivatives along the w direction
     dknots = nurbs.knots;
-    dcoefs = reshape(nurbs.coefs,4*num1*num2,num3);
-    [dcoefs,dknots{3}] = bspderiv(degree(3),dcoefs,nurbs.knots{3});
-    dcoefs = reshape(dcoefs,[4 num1 num2 size(dcoefs,2)]);
-    dnurbs{3} = nrbmak(dcoefs, dknots);
+    dcoefs = reshape (nurbs.coefs,4*num1*num2,num3);
+    [dcoefs,dknots{3}] = bspderiv (degree(3),dcoefs,nurbs.knots{3});
+    dcoefs = reshape (dcoefs,[4 num1 num2 size(dcoefs,2)]);
+    dnurbs{3} = nrbmak (dcoefs, dknots);
 
-  elseif size(nurbs.knots,2) == 2
-  % NURBS structure represents a surface
+    if (nargout == 2)
+      warning ('nrbderiv: the second derivative is not ready for volumes');
+      dnurbs2 = [];
+    end
+
+  elseif (size(nurbs.knots,2) == 2)
+% NURBS structure represents a surface
 
     num1 = nurbs.number(1);
     num2 = nurbs.number(2);
 
-  % taking derivatives along the u direction
+% taking first derivative 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);
+    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);
+
+    if (nargout == 2)
+% taking second derivative along the u direction (duu)
+      dknots2 = dknots;
+      dcoefs2 = permute (dcoefs, [1 3 2]);
+      dcoefs2 = reshape (dcoefs2, 4*num2, []);
+      [dcoefs2, dknots2{1}] = bspderiv (degree(1)-1, dcoefs2, dknots{1});
+      dcoefs2 = permute (reshape (dcoefs2, 4, num2, []), [1 3 2]);
+      dnurbs2{1,1} = nrbmak (dcoefs2, dknots2); 
 
-  % taking derivatives along the v direction
+% taking second derivative along the v direction (duv and dvu)
+      dknots2 = dknots;
+      dcoefs2 = reshape (dcoefs, 4*(num1-1), num2);
+      [dcoefs2, dknots2{2}] = bspderiv (degree(2), dcoefs2, dknots{2});
+      dcoefs2 = reshape (dcoefs2, 4, num1-1, []);
+      dnurbs2{1,2} = nrbmak (dcoefs2, dknots2);
+      dnurbs2{2,1} = dnurbs2{1,2};
+    end
+
+% taking first derivative 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);
+    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);
+
+    if (nargout == 2)
+% taking second derivative along the v direction (dvv)
+      dknots2 = dknots;
+      dcoefs2 = reshape (dcoefs, 4*num1, num2-1);
+      [dcoefs2, dknots2{2}] = bspderiv (degree(2)-1, dcoefs2, dknots{2});
+      dcoefs2 = reshape (dcoefs2, 4, num1, []);
+      dnurbs2{2,2} = nrbmak (dcoefs2, dknots2);
+    end
+
   end
 else
   % NURBS structure represents a curve
 
-  [dcoefs,dknots] = bspderiv(degree,nurbs.coefs,nurbs.knots);
-  dnurbs = nrbmak(dcoefs, dknots);
+  [dcoefs,dknots] = bspderiv (degree, nurbs.coefs, nurbs.knots);
+  dnurbs = nrbmak (dcoefs, dknots);
+  if (nargout == 2)
+    [dcoefs2,dknots2] = bspderiv (degree-1, dcoefs, dknots);
+    dnurbs2 = nrbmak (dcoefs2, dknots2);
+  end
 
 end
 
+varargout{1} = dnurbs;
+if (nargout == 2)
+  varargout{2} = dnurbs2;
+end
+
 end
 
 %!demo

--- a/extra/nurbs/inst/nrbdeval.m	Fri Jul 22 13:22:03 2011 +0000
+++ b/extra/nurbs/inst/nrbdeval.m	Fri Jul 22 15:17:10 2011 +0000
@@ -1,25 +1,21 @@
-function [pnt,jac] = nrbdeval(nurbs, dnurbs, tt)
+function varargout = nrbdeval (nurbs, dnurbs, varargin)
 
-% NRBDEVAL: Evaluation of the derivative NURBS curve, surface or volume.
+% NRBDEVAL: Evaluation of the derivative and second derivatives of NURBS curve, surface or volume.
 %
-%     [pnt, jac] = nrbdeval(crv, dcrv, tt)
-%     [pnt, jac] = nrbdeval(srf, dsrf, {tu tv})
-%     [pnt, jac] = nrbdeval(vol, dvol, {tu tv tw})
+%     [pnt, jac] = nrbdeval (crv, dcrv, tt)
+%     [pnt, jac] = nrbdeval (srf, dsrf, {tu tv})
+%     [pnt, jac] = nrbdeval (vol, dvol, {tu tv tw})
+%     [pnt, jac, hess] = nrbdeval (crv, dcrv, dcrv2, tt)
+%     [pnt, jac, hess] = nrbdeval (srf, dsrf, dsrf2, {tu tv})
+%     [pnt, jac, hess] = nrbdeval (vol, dvol, {tu tv tw})
 %
 % INPUTS:
 %
-%   crv    - original NURBS curve.
-%
-%   srf    - original NUBRS surface
-%
-%   vol    - original NUBRS volume
-%
-%   dcrv   - NURBS derivative representation of crv
-%
-%   dsrf   - NURBS derivative representation of surface
-%
-%   dvol   - NURBS derivative representation of volume
-%
+%   crv,   srf,   vol   - original NURBS curve, surface or volume.
+%   dcrv,  dsrf,  dvol  - NURBS derivative representation of crv, srf 
+%                          or vol (see nrbderiv2)
+%   dcrv2, dsrf2, dvol2 - NURBS second derivative representation of crv,
+%                          srf or vol (see nrbderiv2)
 %   tt     - parametric evaluation points
 %            If the nurbs is a surface or a volume then tt is a cell
 %            {tu, tv} or {tu, tv, tw} are the parametric coordinates
@@ -28,18 +24,11 @@
 %
 %   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);
+%   hess - evaluated second derivatives (Hessian).
 %
 % Copyright (C) 2000 Mark Spink 
-% Copyright (C) 2010 Rafael Vazquez
 % Copyright (C) 2010 Carlo de Falco
+% Copyright (C) 2010, 2011 Rafael Vazquez
 %
 %    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
@@ -54,6 +43,15 @@
 %    You should have received a copy of the GNU General Public License
 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
+if (nargin == 3)
+  tt = varargin{1};
+elseif (nargin == 4)
+  dnurbs2 = varargin{1};
+  tt = varargin{2};
+else 
+  error ('nrbrdeval: wrong number of input parameters')
+end
+
 if (~isstruct(nurbs))
   error('NURBS representation is not structure!');
 end
@@ -62,7 +60,7 @@
   error('Not a recognised NURBS representation');
 end
 
-[cp,cw] = nrbeval (nurbs, tt);
+[cp,cw] = nrbeval(nurbs, tt);
 
 if (iscell(nurbs.knots))
   if (size(nurbs.knots,2) == 3)
@@ -82,6 +80,11 @@
     tempw = cww(ones(3,1),:,:,:);
     jac{3} = (cwp-tempw.*pnt)./temp;
 
+    if (nargout == 3)
+      warning ('nrbdeval: the second derivative is not ready for volumes');
+      hess = [];
+    end
+
   elseif (size(nurbs.knots,2) == 2)
   % NURBS structure represents a surface
     temp = cw(ones(3,1),:,:);
@@ -95,6 +98,22 @@
     tempv = cvw(ones(3,1),:,:);
     jac{2} = (cvp-tempv.*pnt)./temp;
 
+% second derivatives
+    if (nargout == 3 && exist ('dnurbs2'))
+      [cuup, cuuw] = nrbeval (dnurbs2{1,1}, tt);
+      tempuu = cuuw(ones(3,1),:,:);
+      hess{1,1} = (cuup - (2*cup.*tempu + cp.*tempuu)./temp + 2*cp.*tempu.^2./temp.^2)./temp;
+
+      [cvvp, cvvw] = nrbeval (dnurbs2{2,2}, tt);
+      tempvv = cvvw(ones(3,1),:,:);
+      hess{2,2} = (cvvp - (2*cvp.*tempv + cp.*tempvv)./temp + 2*cp.*tempv.^2./temp.^2)./temp;
+
+      [cuvp, cuvw] = nrbeval (dnurbs2{1,2}, tt);
+      tempuv = cuvw(ones(3,1),:,:);
+      hess{1,2} = (cuvp - (cup.*tempv + cvp.*tempu + cp.*tempuv)./temp + 2*cp.*tempu.*tempv./temp.^2)./temp;
+      hess{2,1} = hess{1,2};
+    end
+
   end
 else
 
@@ -107,6 +126,19 @@
   temp1 = cuw(ones(3,1),:);
   jac = (cup-temp1.*pnt)./temp;
 
+  % second derivative
+  if (nargout == 3 && exist ('dnurbs2'))
+    [cuup,cuuw] = nrbeval (dnurbs2, tt);
+    temp2 = cuuw(ones(3,1),:);
+    hess = (cuup - (2*cup.*temp1 + cp.*temp2)./temp + 2*cp.*temp1.^2./temp.^2)./temp;
+  end
+
+end
+
+varargout{1} = pnt;
+varargout{2} = jac;
+if (nargout == 3)
+  varargout{3} = hess;
 end
 
 end