3D B-Spline approximation

Question

B-Spline in 3D can be extended to produce a scalar function of three parameters:

The B-Spline is defined by 64 control points (the data values within a 4X4X4 voxel neighborhood), and evaluated inside the unit cube bounded by the eight central voxels, with t, s and r representing distances in the x, y and z direction respectively.

based on above notes I wrote my function in Matlab:

function [TO]=bspline_matrixform(I,valr,vals,valt,i,j,k)

% I :  Input
% i :  Index of grid point in x direction
% j :  Index of grid point in y direction
% k :  Index of grid point in z direction

r = valr; % Val r value in the grid
s = vals; % Val s value in the grid
t = valt; % Val t value in the grid

M = 1/6*[-1 3 -3 1;3 -6 3 0;-3 0 3 0;1 4 1 0];

R = [ r^3 r^2 r 1]; % 0<=r<=1
S = [ s^3 s^2 s 1]; % 0<=s<=1
T = [ t^3 t^2 t 1]; % 0<=w<=1


index = 1;
QK = zeros(4,1);
for ki= k-1:k+2
    Q = QCalculation(I,i,j,ki);
    QK(index) = S*M*Q*M'*T';
    index = index+1;
end
TO = R*M*[QK(1) QK(2) QK(3) QK(4)]';  

end

function Q = QCalculation (I,i,j,k)

Q = [I(i-1,j-1,k)  I(i,j-1,k)  I(i+1,j-1,k)  I(i+2,j-1,k)
     I(i-1,j,k)    I(i,j,k)    I(i+1,j,k)    I(i+2,j,k)
     I(i-1,j+1,k)  I(i,j+1,k)  I(i+1,j+1,k)  I(i+2,j+1,k)
     I(i-1,j+2,k)  I(i,j+2,k)  I(i+1,j+2,k)  I(i+2,j+2,k)];

end

if my function is true, how should I call my function now? I dont know how should i calculate r,s and t.

P.S: to run this function, i have control points and ngrid in 3D, E.g from X-Y point of view, my control point( blue) and grid(red) are like this:

based on above notes, s,r and t are distances values between control points and the grid in above picture?

Anybody can help me to find out am I going through the right way and understanding? any alternative way or example code that can help im in 3D Bspline?

P.S: I am doing image registration (3d MRI), and my control points are my motion fields of registration result. I need to apply B-Spline to make my registration more smooth.

Thanks a lot