function d2 = dtiXformTensors(d, x)
% dt6_new = dtiXformTensors(dt6, xform)
% 
% Efficiently applies the 3x3 transform to a tensor volume in dt6 format.
% dt6 is a XxYxZx3x3 array.
%
% This is just a matlab-efficient way to do xform*T*xform' on a
% bunch of tensors (T) in dt6 format, a 3d array of the 6 unique
% tensor elements in [Dxx, Dyy, Dzz, Dxy, Dxz, Dyz] order.
% 
%
% 
% HISTORY:
% 2003.12.15 RFD & AHS Wrote it.
% 2006.08.03 RFD: rewote it to be much more memory efficient. It's
% also now ~2.5x faster.

if(size(x,2)>3)
  x = x(1:3,1:3);
end

d2 = zeros(size(d));
d2(:,:,:,1) = x(1)*(x(1).*d(:,:,:,1) + x(4).*d(:,:,:,4) + x(7).*d(:,:,:,5)) ...
            + x(4)*(x(1).*d(:,:,:,4) + x(4).*d(:,:,:,2) + x(7).*d(:,:,:,6)) ...
	    + x(7)*(x(1).*d(:,:,:,5) + x(4).*d(:,:,:,6) + x(7).*d(:,:,:,3));

d2(:,:,:,2) = x(2)*(x(2).*d(:,:,:,1) + x(5).*d(:,:,:,4) + x(8).*d(:,:,:,5)) ...
            + x(5)*(x(2).*d(:,:,:,4) + x(5).*d(:,:,:,2) + x(8).*d(:,:,:,6)) ...
	    + x(8)*(x(2).*d(:,:,:,5) + x(5).*d(:,:,:,6) + x(8).*d(:,:,:,3));

d2(:,:,:,3) = x(3)*(x(3).*d(:,:,:,1) + x(6).*d(:,:,:,4) + x(9).*d(:,:,:,5)) ...
            + x(6)*(x(3).*d(:,:,:,4) + x(6).*d(:,:,:,2) + x(9).*d(:,:,:,6)) ...
	    + x(9)*(x(3).*d(:,:,:,5) + x(6).*d(:,:,:,6) + x(9).*d(:,:,:,3));

d2(:,:,:,4) = x(2)*(x(1).*d(:,:,:,1) + x(4).*d(:,:,:,4) + x(7).*d(:,:,:,5)) ...
            + x(5)*(x(1).*d(:,:,:,4) + x(4).*d(:,:,:,2) + x(7).*d(:,:,:,6)) ...
	    + x(8)*(x(1).*d(:,:,:,5) + x(4).*d(:,:,:,6) + x(7).*d(:,:,:,3));

d2(:,:,:,5) = x(3)*(x(1).*d(:,:,:,1) + x(4).*d(:,:,:,4) + x(7).*d(:,:,:,5)) ...
            + x(6)*(x(1).*d(:,:,:,4) + x(4).*d(:,:,:,2) + x(7).*d(:,:,:,6)) ...
	    + x(9)*(x(1).*d(:,:,:,5) + x(4).*d(:,:,:,6) + x(7).*d(:,:,:,3));

d2(:,:,:,6) = x(3)*(x(2).*d(:,:,:,1) + x(5).*d(:,:,:,4) + x(8).*d(:,:,:,5)) ...
            + x(6)*(x(2).*d(:,:,:,4) + x(5).*d(:,:,:,2) + x(8).*d(:,:,:,6)) ...
	    + x(9)*(x(2).*d(:,:,:,5) + x(5).*d(:,:,:,6) + x(8).*d(:,:,:,3));

return;


% OLD CODE (pre 2006.08.03 rewrite):

temp = zeros([size(d,1), size(d,2), size(d,3), 3, 3]);
temp(:,:,:,1,1) = x(1,1).*d(:,:,:,1) + x(1,2).*d(:,:,:,4) + x(1,3).*d(:,:,:,5);
temp(:,:,:,1,2) = x(1,1).*d(:,:,:,4) + x(1,2).*d(:,:,:,2) + x(1,3).*d(:,:,:,6);
temp(:,:,:,1,3) = x(1,1).*d(:,:,:,5) + x(1,2).*d(:,:,:,6) + x(1,3).*d(:,:,:,3);
temp(:,:,:,2,1) = x(2,1).*d(:,:,:,1) + x(2,2).*d(:,:,:,4) + x(2,3).*d(:,:,:,5);
temp(:,:,:,2,2) = x(2,1).*d(:,:,:,4) + x(2,2).*d(:,:,:,2) + x(2,3).*d(:,:,:,6);
temp(:,:,:,2,3) = x(2,1).*d(:,:,:,5) + x(2,2).*d(:,:,:,6) + x(2,3).*d(:,:,:,3);
temp(:,:,:,3,1) = x(3,1).*d(:,:,:,1) + x(3,2).*d(:,:,:,4) + x(3,3).*d(:,:,:,5);
temp(:,:,:,3,2) = x(3,1).*d(:,:,:,4) + x(3,2).*d(:,:,:,2) + x(3,3).*d(:,:,:,6);
temp(:,:,:,3,3) = x(3,1).*d(:,:,:,5) + x(3,2).*d(:,:,:,6) + x(3,3).*d(:,:,:,3);

temp = dtiXformVectors(temp, x', 'post');
d2 = temp(:,:,:,[1 5 9 4 7 8]);

% This is equivalent to:
% dt6_new = zeros(size(dt6));
% for (x=1:size(dt6,1)),
%   for (y=1:size(dt6,2)),
%       for (z=1:size(dt6,3)),
%           t = dti6to33(dt6_new(x,y,z,:));
%           t = xform * t * xform';
%           dt6_new(x,y,z,:) = dti33to6(t);
%       end
%   end
% end
