Code covered by the BSD License

# Spline2D or Piecewise Continuous 2D Polynomials

### Mark Mikofski (view profile)

Fit a 2D function with piecewise continuous polynomials

splineFit2D(f,x,y,n,m,xb,yb,continuity,ngrid)
```function pp = splineFit2D(f,x,y,n,m,xb,yb,continuity,ngrid)
%SPLINEFIT2D Fit 2D data to a set of piecewise continuous polynomials.
%   PP = SPLINEFIT2D(F,X,Y,N,M,XB,YB) fits F(X,Y) with polynomials of order
%   N, M respectively, within the each pair of consecutive bounds [XB(i) XB(i+1)]
%   and [YB(i) YB(i + 1)] for i from 1 to NUMEL(XB)-1 and NUMEL(YB)-1, respectively.

%% check inputs args
if nargin<8
continuity = 0; % C0 is default continuity
end
gridFlag = false;
if nargin<9 || isempty(ngrid)
gridFlag = true; % use grid spacing determined by number of points in sliver
end
% TODO: check inputs
%% initiallization
% convert all inputs to linear indexing (column vectors)
x = x(:);y = y(:);f = f(:);
Npoints = numel(f); % number of data points
assert(all([numel(x),numel(y)]==Npoints),'splineFit2D:sizeMismatch', ...
'X,Y and F must be same size.')
% number of breaks
Nxb = numel(xb); % x-breaks
Nyb = numel(yb); % y-breaks
xmin = min(x);xmax = max(x);ymin = min(y);ymax = max(y);
xb = [xmin;xb(:);xmax]; % x-breaks
yb = [ymin;yb(:);ymax]; % y-breaks
assert(all([[size(n,1),size(m,1)]==Nxb+1,[size(n,2),size(m,2)]==Nyb+1]), ...
'splineFit2D:sizeMismatch','N and M must have size [numel(xb),numel(yb)].')
% We use data to store all data initially. As we iterated through the
% regions, we remove those points from data that were in that sliver, so
% that only the unsliced data is left - until all of the points have been
% sliced into a sliver separated by breaks.
data = [x y]; % array of all data
% We need to sort F(X,Y) according the breaks, so we use the same method as
% what we're doing for data.
fdata = f; % store F(X,Y) for slicing and sorting
% We store all of the values in cells because we are going to create a
% sparse matrix with all of the data, but we don't know how many data are
% in each sliver until we slice it.
fz = cell(Nxb+1,Nyb+1); % F sorted according to breaks
fi = cell(Nxb+1,Nyb+1); % indices of F in FZ
fidx = 1; % index pf FZ in F, used to sort F according to breaks
z = cell(Nxb+1,Nyb+1); % matrices of terms separated according to breaks
% z indices [i j] in sparse matrix
zi = cell(Nxb+1,Nyb+1); % row
zj = cell(Nxb+1,Nyb+1); % column
zidx = 1;
% matrices of continuity terms at break boundaries
zx0 = cell(Nxb,Nyb+1); % x-breaks
zy0 = cell(Nxb+1,Nyb); % y-breaks
zx0i = cell(Nxb,Nyb+1); % x-breaks
zx0j = cell(Nxb,Nyb+1); % x-breaks
z0idx = 1+Npoints;
xCoeffIdx = 1+sum((n(1,:)+1).*(m(1,:)+1)); % number of coefficients in 1st slice
zy0i = cell(Nxb+1,Nyb); % y-breaks
zy0j = cell(Nxb+1,Nyb); % y-breaks
% higher order continuity requirements
switch continuity
case 0
% this is the default
case 1
%         z1 = cell(Nxb+1,Nyb+1);
case 2
%         z2 = cell(Nxb+1,Nyb+1);
otherwise
error('poop')
end
%% loop over breaks
for xbi = 1:Nxb % x-breaks
% this loop slices the 2D domain into "slices" in the x-direction
%% get indices of slice data
% x-boundaries of slice
xlb = xb(xbi); % lower
xub = xb(xbi+1); % upper
xbIdx = data(:,1)<xub; % indices of x below x-break point
dataSlice = data(xbIdx,:); % slice of data below x-break point
fslice = fdata(xbIdx); % corresponding F(X,Y) of slice
%% loop over y-breaks
for ybi = 1:Nyb % y-breaks
% this loop calculates terms and continuity conditions in each
% "sliver" in the y-direction of the the slice
%% get indices of sliver data
% y-boundaries of sliver
ylb = yb(ybi); % lower
yub = yb(ybi+1); % upper
ybIdx = dataSlice(:,2)<yub; % indices of sliver of data slice below both x- & y-break points
dataSliver = dataSlice(ybIdx,:); % sliver of data slice below below both x- & y-break points
%% matrix of terms for sliver of data slice
z{xbi,ybi} = hornerLoop(dataSliver,n(xbi,ybi),m(xbi,ybi));
%% matrix of terms between slices (x-direction)
Nsliver = size(dataSliver,1); % number of points in sliver
if gridFlag
ngrid = Nsliver; % user did not specify grid spacing at boundaries
end
gridData = [xub*ones(ngrid,1) linspace(ylb,yub,ngrid)']; % grid betwen slices
zx0{xbi,ybi} = [hornerLoop(gridData,n(xbi,ybi),m(xbi,ybi)), ... % left side of slice equals
-hornerLoop(gridData,n(xbi+1,ybi),m(xbi+1,ybi))]; % right side of slice
%% matrix of terms between slivers (y-direction)
gridData = [linspace(xlb,xub,ngrid)' yub*ones(ngrid,1)]; % grid betwen slivers
zy0{xbi,ybi} = [hornerLoop(gridData,n(xbi,ybi),m(xbi,ybi)), ... % bottom side of sliver equals
-hornerLoop(gridData,n(xbi,ybi+1),m(xbi,ybi+1))]; % top side of sliver
dataSlice(ybIdx,:) = []; % remove sliver from slice
%% store F(X,Y) for sliver
fz{xbi,ybi} = fslice(ybIdx);
fi{xbi,ybi} = fidx:fidx+Nsliver-1; % indices
fidx = fidx+Nsliver; % adjust index to next sliver
fslice(ybIdx) = []; % remove sliver of F(X,Y) from slice
%% indices
Ncoeffs = (n(xbi,ybi)+1)*(m(xbi,ybi)+1);
zi{xbi,ybi} = fi{xbi,ybi}'*ones(1,Ncoeffs);
zj{xbi,ybi} = ones(Nsliver,1)*(zidx:zidx+Ncoeffs-1);
zidx = zidx+Ncoeffs;
%% indices of x-boundary
Nxcoeffs = (n(xbi+1,ybi)+1)*(m(xbi+1,ybi)+1);
zx0i{xbi,ybi} = (z0idx:z0idx+ngrid-1)'*ones(1,Ncoeffs+Nxcoeffs);
z0idx = z0idx+ngrid;
zx0j{xbi,ybi} = [zj{xbi,ybi} ones(ngrid,1)*(xCoeffIdx:xCoeffIdx+Nxcoeffs-1)];
xCoeffIdx = xCoeffIdx+Nxcoeffs;
%% indices of y-boundary
Nycoeffs = (n(xbi,ybi+1)+1)*(m(xbi,ybi+1)+1);
zy0i{xbi,ybi} = (z0idx:z0idx+ngrid-1)'*ones(1,Ncoeffs+Nycoeffs);
z0idx = z0idx+ngrid;
zy0j{xbi,ybi} = [zj{xbi,ybi} ones(ngrid,1)*(zidx:zidx+Nycoeffs-1)];
end
%% matrix of terms for last sliver in data slice
% y-boundaries of sliver
ylb = yb(Nyb+1); % lower
yub = yb(Nyb+2); % upper
% after iterating over all of the slivers in the slice, there is still
% data above the last y-break. the remaining  sliver is DATASLICE because as
% we iterated, we removed each sliver from the slice.
z{xbi,Nyb+1} = hornerLoop(dataSlice,n(xbi,Nyb+1),m(xbi,Nyb+1));
%% matrix of terms between slices (x-direction)
% in the last sliver there is only a boundary between neighboring
% slices in the x-direction.
Nsliver = size(dataSlice,1); % number of points in last sliver of slice
if gridFlag
ngrid = Nsliver; % user did not specify grid spacing at boundaries
end
gridData = [xub*ones(ngrid,1) linspace(ylb,yub,ngrid)']; % grid betwen slices
zx0{xbi,Nyb+1} = [hornerLoop(gridData,n(xbi,Nyb+1),m(xbi,Nyb+1)), ...
-hornerLoop(gridData,n(xbi+1,Nyb+1),m(xbi+1,Nyb+1))];
data(xbIdx,:) = []; % remove data slice
%% store F(X,Y) for last sliver of slice
fz{xbi,Nyb+1} = fslice;
fi{xbi,Nyb+1} = fidx:fidx+Nsliver-1; % indices
fidx = fidx+Nsliver; % adjust index to next sliver
fdata(xbIdx) = []; % remove last slice of F(X,Y) from data
%% indices
Ncoeffs = (n(xbi,Nyb+1)+1)*(m(xbi,Nyb+1)+1);
zi{xbi,Nyb+1} = fi{xbi,Nyb+1}'*ones(1,Ncoeffs);
zj{xbi,Nyb+1} = ones(Nsliver,1)*(zidx:zidx+Ncoeffs-1);
zidx = zidx+Ncoeffs;
%% indices of x-boundary
Nxcoeffs = (n(xbi+1,Nyb+1)+1)*(m(xbi+1,Nyb+1)+1);
zx0i{xbi,Nyb+1} = (z0idx:z0idx+ngrid-1)'*ones(1,Ncoeffs+Nxcoeffs);
z0idx = z0idx+ngrid;
zx0j{xbi,Nyb+1} = [zj{xbi,Nyb+1} ones(ngrid,1)*(xCoeffIdx:xCoeffIdx+Nxcoeffs-1)];
xCoeffIdx = xCoeffIdx+Nxcoeffs;
end
%% last slice
% x-boundaries
xlb = xb(Nxb+1); % lower
xub = xb(Nxb+2); % upper
% the last slice is whatever is left in data, because as we iterated we
% removed each slice from data
dataSlice = data; % remaining slice
fslice = fdata; % remaining F(X,Y)
%% loop over y-breaks in last slice
for ybi = 1:Nyb
%% get indices of sliver data
yub = yb(ybi+1); % upper y-bound of sliver
ybIdx = dataSlice(:,2)<yub; % indices of sliver data slice below both x- & y-break points
dataSliver = dataSlice(ybIdx,:); % sliver of data slice below below both x- & y-break points
z{Nxb+1,ybi} = hornerLoop(dataSliver,n(Nxb+1,ybi),m(Nxb+1,ybi)); % matrix of terms for sliver of data slice
%% matrix of terms between slivers (y-direction)
% in the last slice there are only continuity conditions between
% slivers in y-direction
Nsliver = size(dataSliver,1); % number of points in sliver
if gridFlag
ngrid = Nsliver; % user did not specify grid spacing at boundaries
end
gridData = [linspace(xlb,xub,ngrid)' yub*ones(ngrid,1)]; % grid betwen slivers
zy0{Nxb+1,ybi} = [hornerLoop(gridData,n(xbi,ybi),m(xbi,ybi)), ...
-hornerLoop(gridData,n(xbi,ybi+1),m(xbi,ybi+1))];
dataSlice(ybIdx,:) = []; % remove sliver
%% store F(X,Y) for sliver in last slice
fz{Nxb+1,ybi} = fslice(ybIdx);
fi{Nxb+1,ybi} = fidx:fidx+Nsliver-1; % indices
fidx = fidx+Nsliver; % adjust index to next sliver
fslice(ybIdx) = []; % remove slice of F(X,Y) from data
%% indices
Ncoeffs = (n(Nxb+1,ybi)+1)*(m(Nxb+1,ybi)+1);
zi{Nxb+1,ybi} = fi{Nxb+1,ybi}'*ones(1,Ncoeffs);
zj{Nxb+1,ybi} = ones(Nsliver,1)*(zidx:zidx+Ncoeffs-1);
zidx = zidx+Ncoeffs;
%% indices of y-boundary
Nycoeffs = (n(Nxb+1,ybi+1)+1)*(m(Nxb+1,ybi+1)+1);
zy0i{Nxb+1,ybi} = (z0idx:z0idx+ngrid-1)'*ones(1,Ncoeffs+Nycoeffs);
z0idx = z0idx+ngrid;
zy0j{Nxb+1,ybi} = [zj{Nxb+1,ybi} ones(ngrid,1)*(zidx:zidx+Nycoeffs-1)];
end
%% last sliver of data
z{Nxb+1,Nyb+1} = hornerLoop(dataSlice,n(Nxb+1,Nyb+1),m(Nxb+1,Nyb+1)); % matrix of terms for sliver of data slice
% there are no boundaries to match in the last sliver of the last slice
%% store F(X,Y) for sliver in last slice
fz{Nxb+1,Nyb+1} = fslice;
fi{Nxb+1,Nyb+1} = fidx:Npoints; % indices
%% indices
Ncoeffs = (n(Nxb+1,Nyb+1)+1)*(m(Nxb+1,Nyb+1)+1);
zi{Nxb+1,Nyb+1} = fi{Nxb+1,Nyb+1}'*ones(1,Ncoeffs);
zj{Nxb+1,Nyb+1} = ones(Npoints-fidx+1,1)*(zidx:zidx+Ncoeffs-1);
zidx = zidx+Ncoeffs-1; % total number of coefficients
%% fill in matrix of terms for all breaks in sparse matrix
ii = [cellmat2ind(zi);cellmat2ind(zx0i);cellmat2ind(zy0i)];
jj = [cellmat2ind(zj);cellmat2ind(zx0j);cellmat2ind(zy0j)];
s = [cellmat2ind(z);cellmat2ind(zx0);cellmat2ind(zy0)];
mm = z0idx;
nn = zidx;
nzmax = numel(s);
Z = sparse(ii,jj,s,mm,nn,nzmax);
fz = fz';
F = [vertcat(fz{:});zeros(z0idx-Npoints,1)];
form = (n+1).*(m+1);
pp = reshape(mat2cell(Z\F,form(:),1),Nxb+1,Nyb+1);
end

function z = hornerLoop(data,n,m)
%Z = HORNERLOOP(DATA,N,M)
Npoints = size(data,1);
x = data(:,1)*ones(1,n+1);
y = data(:,2)*ones(1,(n+1)*(m+1));
z = ones(Npoints,(n+1)*(m+1));
for ni = 1:n
z(:,1:ni) = z(:,1:ni).*x(:,1:ni);
end
for mi = 1:m
mj = (n+1)*mi;
z(:,1:mj) = z(:,1:mj).*y(:,1:mj);
for ni = 1:n
z(:,mj+1:mj+ni) = z(:,mj+1:mj+ni).*x(:,1:ni);
end
end
end

function ind = cellmat2ind(c)
c = cellfun(@(x)x(:),c,'UniformOutput',false);
ind = vertcat(c{:});
end```