Big_Transform(V,p)

function [Vi,M] = Big_Transform(V,p)
%  Big_Transform: rigid body transformation of a big dataset
%
% [Vi,M] = Big_Transform(V,p)
%
% V: Volume to transform, can be integer uint8.
% p: [ti tj tk ai aj ak si sj sk];
%   (translation, rotation, scaling)
% Vi: new interpolated volume using trilinear interpolation
% M: transformation matrix in homogeneous coordinates
% if only one output variable: returns the transformation matrix only
% 
%
% Olivier Salvado, Case Western Reserve University, 20-jun-06

%% Size of teh interpolated volume
%
V = padarray(V,[1 1 1],'replicate','both');
[nr,nc,ns] = size(V);

%% compute matrices
%
T = eye(4);
T(1,4) = -p(1);
T(2,4) = -p(2);
T(3,4) = -p(3);
R1 = eye(4);
R2 = R1;
R3 = R1;
p(4:6) = p(4:6) *pi/180;
R1(2,2) = cos(p(4));
R1(2,3) = sin(p(4));
R1(3,2) = -sin(p(4));
R1(3,3) = cos(p(4));

R2(1,1) = cos(p(5));
R2(3,1) = -sin(p(5));
R2(1,3) = sin(p(5));
R2(3,3) = cos(p(5));

R3(1,1) = cos(p(6));
R3(1,2) = sin(p(6));
R3(2,1) = -sin(p(6));
R3(2,2) = cos(p(6));

S = eye(4);
if length(p) == 9, % scaling is used
    S(1,1) = p(7);
    S(2,2) = p(8);
    S(3,3) = p(9);
end

% --- combine the matrices
M = T*S*R1*R2*R3;

%% check if only matrix is requested
if nargin==1,
    Vi = M;
    return
end

%% go through all the voxels
%
Vi = 0*V;
Xcenter = [nr nc ns -1]/2+0.5;
for r = 2:nr-1,
    for c= 2:nc-1,
        for s=2:ns-1,

            % --- where to interpolate
            X = M*([r c s 1]-Xcenter)'+Xcenter';

            % --- find scaling for linear interpolation
            Xf = floor(X);
            Xc = Xf+1;
            ratio = (X-Xf);

            % --- actual trilinear interpolation
            if Xf(1)<1 || Xc(1)>nr || Xf(2)<1 || Xc(2)>nc || Xf(3)<1 || Xc(3)>ns,
                Vi(r,c,s) = nan;
            else

                Vi(r,c,s) = ...
                    V(Xf(1),Xf(2),Xf(3))*(1-ratio(1))*(1-ratio(2))*(1-ratio(3))+...
                    V(Xc(1),Xf(2),Xf(3))*(ratio(1))*(1-ratio(2))*(1-ratio(3))+...
                    V(Xf(1),Xc(2),Xf(3))*(1-ratio(1))*(ratio(2))*(1-ratio(3))+...
                    V(Xf(1),Xf(2),Xc(3))*(1-ratio(1))*(1-ratio(2))*(ratio(3))+...
                    V(Xc(1),Xf(2),Xc(3))*(ratio(1))*(1-ratio(2))*(ratio(3))+...
                    V(Xf(1),Xc(2),Xc(3))*(1-ratio(1))*(ratio(2))*(ratio(3))+...
                    V(Xc(1),Xc(2),Xf(3))*(ratio(1))*(ratio(2))*(1-ratio(3))+...
                    V(Xc(1),Xc(2),Xc(3))*(ratio(1))*(ratio(2))*(ratio(3));
            end

        end
    end
end


%% remove padding
%
Vi = Vi(2:end-1,2:end-1,2:end-1);