#! /bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh 'READ.ME' <<'END_OF_FILE'
X ***************************************************************************
X * All the software contained in this library is protected by copyright. *
X * Permission to use, copy, modify, and distribute this software for any *
X * purpose without fee is hereby granted, provided that this entire notice *
X * is included in all copies of any software which is or includes a copy *
X * or modification of this software and in all copies of the supporting *
X * documentation for such software. *
X ***************************************************************************
X * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED *
X * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY *
X * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", *
X * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY *
X * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF *
X * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. *
X ***************************************************************************
X * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE *
X * ABOVE STATEMENT. *
X ***************************************************************************
X
X AUTHOR:
X
X LEIF KOBBELT
X UNIVERSITY OF ERLAGEN-NURNBERG, GERMANY
X E-MAIL: Leif.Kobbelt@informatik.uni-erlangen.de
X
X REFERENCE:
X
X - STABLE EVALUATION OF BOX-SPLINES
X NUMERICAL ALGORITHMS, 14 (1997), PP. 377-382
X
X SOFTWARE REVISION DATE:
X
X NOVEMBER 7, 1996
X
X SOFTWARE LANGUAGE:
X
X MATLAB
X
END_OF_FILE
if test 1743 -ne `wc -c <'READ.ME'`; then
echo shar: \"'READ.ME'\" unpacked with wrong size!
fi
# end of 'READ.ME'
fi
if test -f 'box_demo.m' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'box_demo.m'\"
else
echo shar: Extracting \"'box_demo.m'\" \(2056 characters\)
sed "s/^X//" >'box_demo.m' <<'END_OF_FILE'
Xecho off
X%%
X%% Demo-script for the Box-Spline Evaluation (Leif Kobbelt - 11/07/96)
X%%
X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
Xclc, home
X
Xecho on
X
X%% This is a script file to demonstrate the use of the box_eval() function
X%% which allows the stable and efficient evaluation of box-splines.
X
X%% A box-spline is a piecewise polynomial function from R^s to R.
X%% It is defined by a matrix of k mutually distinct (s-dim) column vectors.
X%% The first example uses s = k = 2 ...
X
XX = [ 1 0 ;
X 0 1 ];
X
X%% The multiple occurrence of identical column is logged in a vector nu
X%% holding the number of times each direction is repeated.
X
Xnu = [3;3];
X
Xpause %% press any key to continue
X
X%% The box-spline evaluation can be applied to a whole set of sample
X%% points. Here we construct a uniform grid over an appropriate interval.
X
X[xx,yy] = meshgrid(((1:20)-2)/6,((1:20)-2)/6);
X p = [xx(:) yy(:)];
X
X%% To also demonstrate the efficiency of the algorithm, we measure the
X%% execution time.
X
Xpause %% press any key to continue
X
Xtic
Xb = box_eval(X,nu,p);
Xtoc
X
Xsurf(reshape(b,20,20))
X
X%% The surface you see right now is the graph of a biquadratic tensor
X%% product Box-spline function.
X
Xpause %% press any key to see more examples
X
Xclc, home
X
X%% The next example has k = 4 distinct directions in R^s = R^2
X%% each occurring only once.
X
X X = [ 1 0 1 1 ;
X 0 1 -1 1 ];
X nu = [1;1;1;1];
X[xx,yy] = meshgrid(((1:20)-2)/6,((1:20)-8)/6);
X p = [xx(:) yy(:)];
X
Xtic
Xb = box_eval(X,nu,p);
Xtoc
X
Xsurf(reshape(b,20,20))
X
X%% The resulting function is called the Zwart-Powell-element
X
Xpause %% press any key for one more example
X
Xclc, home
X
X%% The last example uses k = 3 directions in R^s = R^2. Each has
X%% the multiplicity two.
X
X X = [ 1 0 1 ;
X 0 1 1 ];
X nu = [2;2;2];
X[xx,yy] = meshgrid(((1:20)-2)/5,((1:20)-2)/5);
X p = [xx(:) yy(:)];
X
Xtic
Xb = box_eval(X,nu,p);
Xtoc
X
Xsurf(reshape(b,20,20))
X
X%% This last function is obtained by repeating the directions defining
X%% the Courant-element twice.
X
END_OF_FILE
if test 2056 -ne `wc -c <'box_demo.m'`; then
echo shar: \"'box_demo.m'\" unpacked with wrong size!
fi
# end of 'box_demo.m'
fi
if test -f 'box_eval.m' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'box_eval.m'\"
else
echo shar: Extracting \"'box_eval.m'\" \(1781 characters\)
sed "s/^X//" >'box_eval.m' <<'END_OF_FILE'
Xfunction b = box_eval(X,nu,p)
X%%
X%% Efficient and stable evaluation of box-splines (Leif Kobbelt - 5/15/96
X%% - 11/07/96)
X%%
X%% [Vectorized version: Spline is evaluated
X%% at all locations in p simultaneously]
X%%
X%% b = box_eval(X,nu,p)
X%%
X%% b : Returned vector of function values at locations p
X%% X : Matrix of distinct directions defining the box-spline
X%% (columns = s-dim direction vectors)
X%% nu : Vector holding the multiplicities of the columns in X
X%% p : Vector of points where the spline is to be evaluated
X%% (rows = s-dim point vectors)
X%%
X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
X global BoxEv_I % Identity matrix ... eye(BoxEv_k)
X global BoxEv_J % ones(length(BoxEv_p),1) --> matrix - BoxEv_J * row_vector
X global BoxEv_k % Number of rows in BoxEv_X
X global BoxEv_N % Hash table of normal vectors
X global BoxEv_s % Dimension of domain space
X global BoxEv_u % Hashing function
X global BoxEv_p % Vector of points where the spline is to be evaluated
X global BoxEv_X % Matrix of distinct directions (rows) defining the box-spline
X
X BoxEv_p = p;
X BoxEv_X = X';
X
X [BoxEv_k,BoxEv_s] = size(BoxEv_X);
X [ n,BoxEv_s] = size(p);
X
X BoxEv_I = eye(BoxEv_k);
X BoxEv_J = ones(n,1);
X
X%% Hashing function
X
X BoxEv_u = 2.^(0:BoxEv_k-1);
X
X%% Solve normal equation for least norm representation of p
X
X Y = BoxEv_X(find(nu),:);
X Y = (Y'*Y)\Y';
X
X%% Compute hash table of normal vectors
X
X BoxEv_N = box_norm(BoxEv_s-1,BoxEv_k,zeros(1,BoxEv_k));
X
X%% Do recursion ...
X
X b = box_rec(nu,zeros(1,BoxEv_k),Y,p*Y);
X
X%% Garbage Collection
X
Xclear global BoxEv_I BoxEv_J BoxEv_k BoxEv_N BoxEv_s BoxEv_u BoxEv_p BoxEv_X
END_OF_FILE
if test 1781 -ne `wc -c <'box_eval.m'`; then
echo shar: \"'box_eval.m'\" unpacked with wrong size!
fi
# end of 'box_eval.m'
fi
if test -f 'box_norm.m' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'box_norm.m'\"
else
echo shar: Extracting \"'box_norm.m'\" \(1113 characters\)
sed "s/^X//" >'box_norm.m' <<'END_OF_FILE'
Xfunction N = box_norm(t,k,M)
X%%
X%% Precomputation of normal vectors (Leif Kobbelt- 5/15/96
X%% 11/07/96)
X%%
X%% called by box_eval() ... N = box_norm(t,k,M)
X%%
X%% N : Returned hash table of normal vectors
X%% t : Number of rows to be selected before base case is reached
X%% k : Next row of BoxEv_X to be considered for selection
X%% M : Bitvector indicating the selected rows of BoxEv_X
X%%
X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
X global BoxEv_I % Identity matrix ... eye(BoxEv_k)
X global BoxEv_s % Dimension of domain space
X global BoxEv_X % Matrix of distinct directions (rows) defining the box-spline
X
X%% Allocate space according to hashing function
X
X N = zeros(BoxEv_s,2^k);
X
X%% Cut leafless branches
X
X if (k>=t)
X if (t>0)
X
X%% Left half of the hash table: row X(k,:) not selected
X%% Right half of the hash table: row X(k,:) selected
X
X N = [box_norm(t,k-1,M),box_norm(t-1,k-1,M+BoxEv_I(k,:))];
X
X else
X
X%% Normal vector is orthogonal to all selected rows ...
X
X N(:,1) = null(BoxEv_X(find(M),:));
X end
X end
END_OF_FILE
if test 1113 -ne `wc -c <'box_norm.m'`; then
echo shar: \"'box_norm.m'\" unpacked with wrong size!
fi
# end of 'box_norm.m'
fi
if test -f 'box_rec.m' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'box_rec.m'\"
else
echo shar: Extracting \"'box_rec.m'\" \(2625 characters\)
sed "s/^X//" >'box_rec.m' <<'END_OF_FILE'
Xfunction b = box_rec(n,m,Y,t)
X%%
X%% Recursive Evaluation of Box-Splines (Leif Kobbelt - 5/15/96
X%% - 11/07/96)
X%%
X%% called by box_eval() ... b = box_rec(n,m,Y,t)
X%%
X%% b : Returned vector of function values at locations BoxEv_p
X%% n : Vector holding the multiplicities of the rows in BoxEv_X
X%% m : Current position in the recursion tree (for delayed translation)
X%% Y : Matrix to compute the least norm representation of BoxEv_p
X%% with respect to the remaining directions BoxEv_X(find(n),:)
X%% t : Least norm representation of BoxEv_p
X%%
X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X
X global BoxEv_I % Identity matrix ... eye(BoxEv_k)
X global BoxEv_J % ones(length(BoxEv_p),1) --> matrix - BoxEv_J * row_vector
X global BoxEv_k % Number of rows in BoxEv_X
X global BoxEv_N % Hash table of normal vectors
X global BoxEv_s % Dimension of domain space
X global BoxEv_u % Hashing function
X global BoxEv_p % Vector of points where the spline is to be evaluated
X global BoxEv_X % Matrix of distinct directions (rows) defining the box-spline
X
X if (sum(n)>BoxEv_s)
X
X%% Recursion case ...
X
X b = 0;
X j = 1;
X
X%% Sum over the remaining directions in BoxEv_X ...
X
X for i = 1:BoxEv_k
X
X%% Update multiplicity of directions and position in recursion tree
X
X nn = n-BoxEv_I(:,i);
X mm = m+BoxEv_I(i,:);
X
X%% Recursive calls
X
X if (n(i)>1)
X b = b+ t(:,j) .*box_rec(nn,m ,Y,t );
X b = b+(n(i)-t(:,j)).*box_rec(nn,mm,Y,t-BoxEv_J*(BoxEv_X(i,:)*Y));
X j = j+1;
X elseif (n(i)>0)
X
X%% Update least norm representation
X
X Z = BoxEv_X(find(nn),:);
X if (rank(Z) == BoxEv_s)
X Z = (Z'*Z)\Z';
X b = b+ t(:,j) .*box_rec(nn,m ,Z,(BoxEv_p-BoxEv_J*(m *BoxEv_X))*Z);
X b = b+(n(i)-t(:,j)).*box_rec(nn,mm,Z,(BoxEv_p-BoxEv_J*(mm*BoxEv_X))*Z);
X end
X j = j+1;
X end
X end
X
X%% Normalization
X
X b = b/(sum(n)-BoxEv_s);
X else
X
X%% Base case ... compute characteristic function
X
X b = 1;
X
X%% Delayed translations
X
X v = find(n);
X z = BoxEv_p-BoxEv_J*(m*BoxEv_X);
X
X%% Check against all hyperplanes
X
X for i = 1:BoxEv_s
X
X%% Lookup normal vector to current hyperplane
X
X NN = BoxEv_N(:,1+BoxEv_u*(n-BoxEv_I(:,v(i))));
X
X%% Box is half-open!!!
X
X p = BoxEv_X(v(i),:)*NN;
X q = z*NN;
X b = min(b,1-((p>0&q<0)|(p<0&q>=0)));
X q = (BoxEv_p-BoxEv_J*((m+BoxEv_I(v(i),:))*BoxEv_X))*NN;
X b = min(b,1-((p>0&q>=0)|(p<0&q<0)));
X end
X
X%% Normalization
X
X b = b/abs(det(BoxEv_X(v(1:BoxEv_s),:)));
X end
END_OF_FILE
if test 2625 -ne `wc -c <'box_rec.m'`; then
echo shar: \"'box_rec.m'\" unpacked with wrong size!
fi
# end of 'box_rec.m'
fi
echo shar: End of shell archive.
exit 0