Code covered by the BSD License  

Highlights from
Spline2D or Piecewise Continuous 2D Polynomials

image thumbnail

Spline2D or Piecewise Continuous 2D Polynomials

by

 

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

Contact us