function centripetalCatmullRom()

x = [0.1769,0.1009;0.1184,0.2120;0.0716,0.3874;0.0775,0.5629;0.1652,0.6681;0.2792,0.6798;0.3611,0.5629;0.3699,0.4693;0.3699,0.3260;0.4313,0.2120;0.5249,0.1711;0.6477,0.1711;0.7003,0.2354;0.7003,0.3874;0.6740,0.5044;0.6184,0.5716;0.5629,0.6418;0.5629,0.7383;0.6184,0.7675;0.7675,0.7880;0.9225,0.6447;0.9488,0.5278;0.9284,0.4225;0.8494,0.4079;0.7617,0.4284;0.7266,0.4781;0.6915,0.5804;0.6944,0.6243;0.7383,0.6740;0.7880,0.6711]';

% compute centripetal CR-spline intervals
d = x(:, 2:end) - x(:,1:end-1);
dt = sqrt(d(1,:).*d(1,:) + d(2,:).*d(2,:));
tI = [0, cumsum(dt)];
nTI = length(tI);

plot(x(1,:), x(2,:), 'ko');
axis equal;
hold on;

t = linspace(0, tI(end), 400);
itv = compute_intervals(tI, t);

% clip first / last intervals
range = (itv >= 2 & itv < nTI-1);
itv = itv(range);
t = repmat(t(range), [2 1]);

% array setup
t0 = repmat(tI(itv-1), [2 1]);
t1 = repmat(tI(itv), [2 1]);
t2 = repmat(tI(itv+1), [2 1]);
t3 = repmat(tI(itv+2), [2 1]);

% Evaluate CR-Spline
P0 = x(:,itv-1);
P1 = x(:,itv);
P2 = x(:,itv+1);
P3 = x(:,itv+2);

A1 = (t1 - t)./(t1 - t0) .* P0 + (t - t0)./(t1 - t0) .* P1;
A2 = (t2 - t)./(t2 - t1) .* P1 + (t - t1)./(t2 - t1) .* P2;
A3 = (t3 - t)./(t3 - t2) .* P2 + (t - t2)./(t3 - t2) .* P3;

B1 = (t2 - t)./(t2 - t0) .* A1 + (t - t0)./(t2 - t0) .* A2;
B2 = (t3 - t)./(t3 - t1) .* A2 + (t - t1)./(t3 - t1) .* A3;

C = (t2 - t)./(t2 - t1) .* B1 + (t - t1)./(t2 - t1) .* B2;

% Evaluate CR-Spline Derivatives
A1p = (P1 - P0)./(t1 - t0);
A2p = (P2 - P1)./(t2 - t1);
A3p = (P3 - P2)./(t3 - t2);

B1p = (A2 - A1) ./ (t2 - t0) + (t2 - t)./(t2 - t0) .* A1p + (t - t0)./(t2 - t0) .* A2p;
B2p = (A3 - A2) ./ (t3 - t1) + (t3 - t)./(t3 - t1) .* A2p + (t - t1)./(t3 - t1) .* A3p;

Cp = (B2 - B1) ./ (t2 - t1) + (t2 - t)./(t2 - t1) .* B1p + (t - t1)./(t2 - t1) .* B2p;

plot(C(1,:), C(2,:), 'r-');
quiver(C(1,:), C(2,:), Cp(1,:), Cp(2,:), 0.1, 'k-');


% returns interval index for evaluation points
    function iv = compute_intervals(tt, vv)
        idx = 1;
        
        iv = 0*vv;
        
        ivIdx = 1;
        for xx = vv
            while tt(idx+1) < xx
                idx=idx+1;
            end
            
            iv(ivIdx) = idx;
            ivIdx = ivIdx+1;
        end
        
        
    end
end

