Spline2D - tools to fit 2D functions with piecewise continuous polynomials.

examples.m

%% spline2D example
% Initialize a clean workspace.
% Save your work first!
clear all
close all
%% Set up test 2D polynomial.
[x,y] = meshgrid(-3:4,-3:4); % independent variables
% EG: Make a test 2D polynomial with coeffiecients 1, 2, 3, 4, 5, 6
% that is 2nd-order (quadratic) in x-direction, and 1st order (linear)
% in y-direction.
test2DPoly = polyVal2D([1,2,3,4,5,6],x,y,2,1); % a test 2D polynomial
%% Fit test polynomial with spline (piecewise continuous polynomials).
% Split region at 0 in x-direction and at -1 and 1 in y-direction, so there are
% 6 regions and also 6 boundaries to fit. In each of the 6 regions fit
% to 2nd order in x-direction and 1st order in y-direction. Use the default
% number of gridpoints along the boundaries (IE: don't specify).
fit2DSpline = splineFit2D(test2DPoly,x,y,2*ones(2,3),ones(2,3),0,[-1 1]); % spline fit
test2DSpline = splineVal2D(fit2DSpline,x,y,2*ones(2,3),ones(2,3),0,[-1 1]); % evaluate
% another example with 9 regions and 9 boundaries!
% fit2DSpline = splineFit2D(test2DPoly,x,y,2*ones(3),ones(3),[-1 1],[-1 1]);
% test2DSpline = splineVal2D(fit2DSpline,x,y,2*ones(3),ones(3),[-1 1],[-1 1]);
%% plot results
testRange = [min(test2DPoly(:));max(test2DPoly(:))];
testRange = [testRange(1) - diff(testRange)/10; ...
    testRange(2) + diff(testRange)/10] * [1 1];
figure, subplot(1,2,1)
scatter(x(:),test2DPoly(:),750,y(:),'.'), hold('on')
plot(x',test2DSpline')
plot([0 0],testRange,'--','LineWidth',2.5), hold('off'), grid
title('x-direction')
xlabel('x'), ylabel('test2DPoly')
legstr = cellstr([repmat('y_{spline} = ',8,1),num2str((-3:4)')]);
legend('poly',legstr{:})
subplot(1,2,2)
scatter(y(:),test2DPoly(:),750,x(:),'.'), hold('on')
plot(y,test2DSpline)
plot([-1 1;-1 1],testRange,'--','LineWidth',2.5), hold('off'), grid
title('y-direction')
xlabel('y'), ylabel('test2DPoly')
legstr = cellstr([repmat('y_{spline} = ',8,1),num2str((-3:4)')]);
legend('poly',legstr{:})

splineFit2D.m

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
% add extrema to breaks
xmin = min(x);xmax = max(x);ymin = min(y);ymax = max(y);
assert(all(xmin<xb(:)) && all(xb(:)<xmax),'splineFit2D:badBreak','bad break')
assert(all(ymin<yb(:)) && all(yb(:)<ymax),'splineFit2D:badBreak','bad break')
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

splineVal2D.m

function f = splineVal2D(pp,x,y,n,m,xb,yb)
%SPLINEVAL2D Evaluate 2D data to a set of piecewise continuous polynomials.
%   F = SPLINEFIT2D(PP,X,Y,N,M,XB,YB) evaluates 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
% TODO: check inputs
%% initiallization
% convert all inputs to linear indexing (column vectors)
assert(all(size(x)==size(y)),'polyVal2D:sizeMismatch', ...
    'X and Y must be the same size.')
% number of breaks
Nxb = numel(xb); % x-breaks
Nyb = numel(yb); % 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.
f = zeros(size(x)); % F sorted according to breaks
fbool = true(size(x)); % F sorted according to breaks
%% 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
    xbIdx = data(:,1)<xb(xbi); % indices of x below x-break point
    dataSlice = data(xbIdx,:); % slice of data below x-break point
    fslice = zeros(size(dataSlice,1),1);
    fboolslice = true(size(dataSlice,1),1);
    %% 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
        ybIdx = dataSlice(:,2)<yb(ybi); % 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
        fidx = find(fboolslice); % indices of remaining points
        fslice(fidx(ybIdx)) = polyVal2D(pp{xbi,ybi},dataSliver(:,1),dataSliver(:,2),n(xbi,ybi),m(xbi,ybi));
        fboolslice(fidx(ybIdx)) = false; % eliminate sliver
        dataSlice(ybIdx,:) = []; % remove sliver from slice
    end
    %% matrix of terms for last sliver in data slice
    % 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.
    fslice(fboolslice) = polyVal2D(pp{xbi,Nyb+1},dataSlice(:,1),dataSlice(:,2),n(xbi,Nyb+1),m(xbi,Nyb+1));
    fidx = find(fbool);
    f(fidx(xbIdx)) = fslice;
    fbool(fidx(xbIdx)) = false;
    data(xbIdx,:) = []; % remove data slice
end
%% last slice
% the last slice is whatever is left in data, because as we iterated we
% removed each slice from data
dataSlice = data; % remaining slice
fslice = zeros(size(dataSlice,1),1);
fboolslice = true(size(dataSlice,1),1);
%% loop over y-breaks in last slice
for ybi = 1:Nyb
    %% get indices of sliver data
    ybIdx = dataSlice(:,2)<yb(ybi); % 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
    fidx = find(fboolslice); % indices of remaining points
    fslice(fidx(ybIdx)) = polyVal2D(pp{Nxb+1,ybi},dataSliver(:,1),dataSliver(:,2),n(Nxb+1,ybi),m(Nxb+1,ybi)); % matrix of terms for sliver of data slice
    fboolslice(fidx(ybIdx)) = false; % eliminate sliver
    dataSlice(ybIdx,:) = []; % remove sliver
end
%% last sliver of data
fslice(fboolslice) = polyVal2D(pp{Nxb+1,Nyb+1},dataSlice(:,1),dataSlice(:,2),n(Nxb+1,Nyb+1),m(Nxb+1,Nyb+1)); % matrix of terms for sliver of data slice
f(fbool) = fslice;
end