Code covered by the BSD License

# Spline2D or Piecewise Continuous 2D Polynomials

### Mark Mikofski (view profile)

Fit a 2D function with piecewise continuous polynomials

splineVal2D(pp,x,y,n,m,xb,yb)
```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
```