Catmull-Rom cubic spline interpolation [Matlab]

crspline_curve.m

function varargout = crspline_curve(varargin)
%CRSPLINE_CURVE Конструирует двумерную кривую кубическим Catmull-Rom интерполяционным сплайном
%
% Описание:
%   Функция конструирует двумерную кривую по узловым точкам, используя
%   Catmull-Rom интерполяционные кубические сплайны.
%   Данные сплайны могут локально контролироваться и проходят по
%   узловым точкам, а не в выпуклой их оболочке.
%   
%   Литература:
%   [1] Christopher Twigg, Catmull-Rom splines, 2003
%   [2] Catmull, E., and Rom, R. A class of local interpolating splines, 1974
%
%
% Использование:
%   C = crspline_curve(P)
%   C = crspline_curve(P, n)
%   C = crspline_curve(P, tension)
%   C = crspline_curve(P, n, tension)
%   C = crspline_curve(px, py)
%   C = crspline_curve(px, py, n)
%   C = crspline_curve(px, py, tension)
%   C = crspline_curve(px, py, n, tension)
%   [cx, cy] = crspline_curve(...)
%
%
% Входные аргументы:
%   P  -- Матрица X,Y координат опорных точек [PX; PY]
%   px -- Вектор X-координат опорных точек
%   py -- Вектор Y-координат опорных точек
%   n  -- Количество точек кривой
%   tension -- задаёт параметр упругости сплайна [0..1]
%
%
% Выходные аругменты:
%   C  -- Матрица X,Y координат точек кривой
%   cx -- Вектор X-координат точек кривой
%   cy -- Вектор Y-координат точек кривой
%

Inputs = parse_inputs(varargin{:});
validate_inputs(Inputs);

px = [Inputs.px(1) Inputs.px Inputs.px(end)];
py = [Inputs.py(1) Inputs.py Inputs.py(end)];

n = Inputs.n;
tension = Inputs.tension;

cpn = length(px);

% Количество интервалов (отрезков)
ni = cpn - 3;

% Количество точек кривой на каждом из интервалов 1:ni-1
npf = floor(n/ni);

% Количество точек на последнем интервале (то, что осталось)
npl = n - (npf * (ni - 1));

% Количества точек в каждом интервале
nps = [npf*ones(1, ni-1), npl];

% Для удаления по 1-й лишней точке из каждого интервала кроме первого
nps(2:end) = nps(2:end) + 1;

% Здесь храним сплайны для каждого интервала
Cc = cell(1, ni);

% Определяем базисную матрицу для CR-сплайна
Mb = [...
            0,         1,           0,        0; ...
     -tension,         0,     tension,        0; ...
    2*tension, tension-3, 3-2*tension, -tension; ...
     -tension, 2-tension,   tension-2,  tension; ...
    ];

% Номер текущего интервала
ii = 0;
for i = 3:cpn-1
    ii = ii + 1;
    
    % Узловые точки (Pi-2, Pi-1, Pi, Pi+1)
    PX = [px(i-2); px(i-1); px(i); px(i+1)];
    PY = [py(i-2); py(i-1); py(i); py(i+1)];
    
    % Вычисляем сплайн для каждого интервала по заданным узловым точкам
    Cc{ii} = eval_crspline(Mb, PX, PY, nps(ii));
end

% Удаляем каждую первую точку в каждом интервале кроме первого
% Нужно для того, чтобы точки сплайна в контрольных точках  не дублировались
Cc(2:end) = cellfun(@(x) x(:,2:end), Cc(2:end), 'UniformOutput', false);

% Получаем полную сплайн-кривую
C = horzcat([], Cc{:});

if (nargout == 1)
    varargout{1} = C;
elseif (nargout == 2)
    varargout{1} = C(1,:);
    varargout{2} = C(2,:);
end


%==========================================================================
function C = eval_crspline(Mb, PX, PY, np)
%EVAL_CRSPLINE Вычислить значения C-R сплайна для заданного количества точек

u = linspace(0, 1, np);
u2 = u.*u;  % u^2
u3 = u2.*u; % u^3

C = zeros(2, np);

for n = 1:np
    C(:,n) = [1 u(n) u2(n) u3(n)] * Mb * [PX, PY];
end


%==========================================================================
function Inputs = parse_inputs(varargin)
%PARSE_INPUTS

error(nargchk(1, 4, nargin))

n = 101;
tension = 0.5;

switch nargin
    case 1
        P = varargin{1};
    case 2
        if isscalar(varargin{2})
            P = varargin{1};
            
            if (varargin{2} <= 1)
                tension = varargin{2};
            else
                n = varargin{2};
            end
        else
            px = varargin{1};
            py = varargin{2};
        end
    case 3
        if isscalar(varargin{2})
            P = varargin{1};
            n = varargin{2};
            tension = varargin{3};
        else
            px = varargin{1};
            py = varargin{2};
            
            if (varargin{3} <= 1)
                tension = varargin{3};
            else
                n = varargin{3};
            end
        end
    case 4
        px = varargin{1};
        py = varargin{2};
        n = varargin{3};
        tension = varargin{4};
end

if exist('P', 'var')
    if (size(P, 1) ~= 2)
        error('crspline_curve:invalidInputs', ...
            'First input argument must be a matrix of 2 rows.')
    end
    
    px = P(1,:);
    py = P(2,:);
end

Inputs.px = px;
Inputs.py = py;
Inputs.n = n;
Inputs.tension = tension;


%==========================================================================
function validate_inputs(Inputs)
%VALIDATE_INPUTS

validateattributes(Inputs.px, {'numeric'}, {'real', 'nonempty', 'nonnan', 'vector'}, ...
    mfilename('fullpath'), 'PX');
validateattributes(Inputs.py, {'numeric'}, {'real', 'nonempty', 'nonnan', 'vector'}, ...
    mfilename('fullpath'), 'PY');

if (length(Inputs.px) ~= length(Inputs.py))
    error('crspline_curve:invalidInputs', ...
        'Vectors PX and PY must be the same lengths.')
end

validateattributes(Inputs.n, {'numeric'}, {'scalar', '>=', length(Inputs.px)}, ...
    mfilename('fullpath'), 'number of points');

validateattributes(Inputs.tension, {'numeric'}, {'scalar', '>=', 0, '<=', 2}, ...
    mfilename('fullpath'), 'Tension');