#! /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