[Octave-cvsupdate] SF.net SVN: octave:[9785] trunk/octave-forge/extra/nurbs/inst

[<raf...@us...>]

Revision: 9785
          http://octave.svn.sourceforge.net/octave/?rev=9785&view=rev
Author:   rafavzqz
Date:     2012-03-08 16:14:36 +0000 (Thu, 08 Mar 2012)
Log Message:
-----------
Function and derivative computed at the same time, for speedup

Modified Paths:
--------------
    trunk/octave-forge/extra/nurbs/inst/private/onebasisfun__.m
    trunk/octave-forge/extra/nurbs/inst/private/onebasisfunder__.m
    trunk/octave-forge/extra/nurbs/inst/tbasisfun.m

Modified: trunk/octave-forge/extra/nurbs/inst/private/onebasisfun__.m
===================================================================
--- trunk/octave-forge/extra/nurbs/inst/private/onebasisfun__.m	2012-03-08 16:13:11 UTC (rev 9784)
+++ trunk/octave-forge/extra/nurbs/inst/private/onebasisfun__.m	2012-03-08 16:14:36 UTC (rev 9785)
@@ -1,60 +1,56 @@
-function N = onebasisfun__ (u, p, U)
+function Nip = onebasisfun__ (u, p, U)
 
 %  __ONEBASISFUN__: Undocumented internal function
 %
+%   Adapted from Algorithm A2.4 from 'The NURBS BOOK' pg74.
+%
 %   Copyright (C) 2009 Carlo de Falco
 %   Copyright (C) 2012 Rafael Vazquez
 %   This software comes with ABSOLUTELY NO WARRANTY; see the file
 %   COPYING for details.  This is free software, and you are welcome
 %   to distribute it under the conditions laid out in COPYING.
 
-  N = zeros (size (u));
+  Nip = zeros (size (u));
+  N = zeros (p+1, 1);
 
-  for ii = 1:numel (u)
+  for ii = 1:numel(u)
+    if ((u(ii) == U(1)) && (U(1) == U(end-1)) || ...
+        (u(ii) == U(end)) && (U(end) == U(2)))
+      Nip(ii) = 1;
+      continue
+    end
     if (~ any (U <= u(ii))) || (~ any (U > u(ii)))
       continue;
-    elseif (p == 0)
-      N(ii) = 1;
-      continue;
-    elseif (p == 1)
-      if (u(ii) < U(2))
-        N(ii) = (u(ii) - U(1)) / (U(2) - U(1));
-        continue;
+    end
+    for jj = 1:p+1 % Initialize zero-th degree functions
+      if (u(ii) > U(jj) && u(ii) < U(jj+1))
+        N(jj) = 1;
       else
-        N(ii) = (U(end) - u(ii)) / (U(3) - U(2));
-        continue;
+        N(jj) = 0;
       end
-
-    elseif (p == 2)
-      if (u(ii) < U(2))
-        N(ii) = (u(ii) - U(1))^2 / ((U(3) - U(1)) * (U(2) - U(1)));
-        continue;
-      elseif (u(ii) > U(3))
-        N(ii) = (U(4) - u(ii))^2 / ((U(4) - U(2)) * (U(4) - U(3)));
-        continue;
+    end
+    for k = 1:p
+      if (N(1) == 0)
+        saved = 0;
       else
-        ld = U(3) - U(1); dd = U(4) - U(2);
-        if (ld ~= 0)
-          N(ii) = N(ii) + (u(ii) - U(1))*(U(3) - u(ii)) / ((U(3) - U(2)) * ld);
-        end
-        if (dd ~= 0)
-          N(ii) = N(ii) + (U(4) - u(ii))*(u(ii) - U(2)) / ((U(3) - U(2)) * dd);
-        end
+        saved = (u(ii) - U(1))*N(1) / (U(k+1)-U(1));
       end
 
-    else
-      ln = u(ii) - U(1);
-      ld = U(end-1) - U(1);
-      if (ld ~= 0)
-        N(ii) = N(ii) + ln * onebasisfun__ (u(ii), p-1, U(1:end-1))/ ld; 
+      for jj = 1:p-k+1
+        Uleft = U(1+jj);
+        Uright = U(1+jj+k);
+        if (N(jj+1) == 0)
+          N(jj) = saved;
+          saved = 0;
+        else
+          temp = N(jj+1)/(Uright-Uleft);
+          N(jj) = saved + (Uright - u(ii))*temp;
+          saved = (u(ii) - Uleft)*temp;
+        end
       end
-
-      dn = U(end) - u(ii);
-      dd = U(end) - U(2);
-      if (dd ~= 0)
-        N(ii) = N(ii) + dn * onebasisfun__ (u(ii), p-1, U(2:end))/ dd;
-      end
     end
+    Nip(ii) = N(1);
   end
-  
+
+
 end

Modified: trunk/octave-forge/extra/nurbs/inst/private/onebasisfunder__.m
===================================================================
--- trunk/octave-forge/extra/nurbs/inst/private/onebasisfunder__.m	2012-03-08 16:13:11 UTC (rev 9784)
+++ trunk/octave-forge/extra/nurbs/inst/private/onebasisfunder__.m	2012-03-08 16:14:36 UTC (rev 9785)
@@ -1,4 +1,4 @@
-function N = onebasisfunder__ (u, p, U)
+function [N, Nder] = onebasisfunder2__ (u, p, U)
 
 %  __ONEBASISFUNDER__: Undocumented internal function
 %
@@ -8,22 +8,30 @@
 %   to distribute it under the conditions laid out in COPYING.
 
   N = zeros (size (u));
+  Nder = zeros (size (u));
 
   for ii = 1:numel (u)
     if (~ any (U <= u(ii))) || (~ any (U > u(ii)))
       continue;
     elseif (p == 0)
-      N(ii) = 0;
+      N(ii) = 1;
+      Nder(ii) = 0;
       continue;
     else
+      ln = u(ii) - U(1);
       ld = U(end-1) - U(1);
       if (ld ~= 0)
-        N(ii) = N(ii) + p * onebasisfun__ (u(ii), p-1, U(1:end-1))/ ld;
+        aux = onebasisfun__ (u(ii), p-1, U(1:end-1))/ ld;
+        N(ii) = N(ii) + ln * aux;
+        Nder(ii) = Nder(ii) + p * aux;
       end
 
+      dn = U(end) - u(ii);
       dd = U(end) - U(2);
       if (dd ~= 0)
-        N(ii) = N(ii) - p * onebasisfun__ (u(ii), p-1, U(2:end))/ dd;
+        aux = onebasisfun__ (u(ii), p-1, U(2:end))/ dd;
+        N(ii) = N(ii) + dn * aux;
+        Nder(ii) = Nder(ii) - p * aux;
       end
     end
   end

Modified: trunk/octave-forge/extra/nurbs/inst/tbasisfun.m
===================================================================
--- trunk/octave-forge/extra/nurbs/inst/tbasisfun.m	2012-03-08 16:13:11 UTC (rev 9784)
+++ trunk/octave-forge/extra/nurbs/inst/tbasisfun.m	2012-03-08 16:14:36 UTC (rev 9785)
@@ -44,10 +44,10 @@
       error ('tbasisfun: knot vector and degree do not correspond')
     end
     
-    N = onebasisfun__ (u, p, U);
-
-    if (nargout == 2)
-      Nder = onebasisfunder__ (u, p, U);
+    if (nargout == 1)
+      N = onebasisfun__ (u, p, U);
+    else
+      [N, Nder] = onebasisfunder__ (u, p, U);
     end
     
   elseif (size(U,2) == 2)
@@ -56,15 +56,16 @@
       error ('tbasisfun: knot vector and degree do not correspond')
     end
     
-    Nu = onebasisfun__ (u(1,:), p(1), U{1});
-    Nv = onebasisfun__ (u(2,:), p(2), U{2});
+    if (nargout == 1)
+      Nu = onebasisfun__ (u(1,:), p(1), U{1});
+      Nv = onebasisfun__ (u(2,:), p(2), U{2});
 
-    N = Nu.*Nv;
+      N = Nu.*Nv;
+    elseif (nargout == 2)
+      [Nu, Ndu] = onebasisfunder__ (u(1,:), p(1), U{1});
+      [Nv, Ndv] = onebasisfunder__ (u(2,:), p(2), U{2});
 
-    if (nargout == 2)
-      Ndu = onebasisfunder__ (u(1,:), p(1), U{1});
-      Ndv = onebasisfunder__ (u(2,:), p(2), U{2});
-      
+      N = Nu.*Nv;
       Nder(1,:) = Ndu.*Nv;
       Nder(2,:) = Nu.*Ndv;
     end
@@ -75,17 +76,18 @@
       error ('tbasisfun: knot vector and degree do not correspond')
     end
 
-    Nu = onebasisfun__ (u(1,:), p(1), U{1});
-    Nv = onebasisfun__ (u(2,:), p(2), U{2});
-    Nw = onebasisfun__ (u(3,:), p(3), U{3});
+    if (nargout == 1)
+      Nu = onebasisfun__ (u(1,:), p(1), U{1});
+      Nv = onebasisfun__ (u(2,:), p(2), U{2});
+      Nw = onebasisfun__ (u(3,:), p(3), U{3});
 
-    N = Nu.*Nv.*Nw;
-
-    if (nargout == 2)
-      Ndu = onebasisfunder__ (u(1,:), p(1), U{1});
-      Ndv = onebasisfunder__ (u(2,:), p(2), U{2});
-      Ndw = onebasisfunder__ (u(3,:), p(3), U{3});
+      N = Nu.*Nv.*Nw;
+    else
+      [Nu, Ndu] = onebasisfunder__ (u(1,:), p(1), U{1});
+      [Nv, Ndv] = onebasisfunder__ (u(2,:), p(2), U{2});
+      [Nw, Ndw] = onebasisfunder__ (u(3,:), p(3), U{3});
       
+      N = Nu.*Nv.*Nw;
       Nder(1,:) = Ndu.*Nv.*Nw;
       Nder(2,:) = Nu.*Ndv.*Nw;
       Nder(3,:) = Nu.*Nv.*Ndw;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.