function spm_preserve_quantity(matname,infile)
% Modulate spatially normalized images in order to preserve `counts'.
% FORMAT spm_preserve_quantity(matname,infile)
% matname  - the `_sn3d.mat' file containing the spatial normalization
%            parameters.
% infile   - the spatially normalized image to be modulated.
%
% After nonlinear spatial normalization, the relative volumes of some
% brain structures will have decreased, whereas others will increase.
% The resampling of the images preserves the concentration of pixel
% units in the images, so the total counts from structures that have
% reduced volumes after spatial normalization will be reduced by an
% amount proportional to the volume reduction.
%
% This routine rescales images after spatial normalization, so that
% the total counts from any structure are preserved.  It was written
% as an optional step in performing voxel based morphometry.
%
%_______________________________________________________________________
% %W% John Ashburner %E%

load(deblank(matname))
if (exist('mgc') ~= 1)
	error(['Matrix file ' matname ' is the wrong type.']);
end
if (mgc ~= 960209)
	error(['Matrix file ' matname ' is the wrong type.']);
end

V      = spm_vol(infile);

% Assume transverse images, and obtain
% position of pixels in millimeters.
%----------------------------------------------------------------------------
x = (1:V.dim(1))*V.mat(1,1) + V.mat(1,4);
y = (1:V.dim(2))*V.mat(2,2) + V.mat(2,4);
z = (1:V.dim(3))*V.mat(3,3) + V.mat(3,4);

% Convert to voxel space of template.
%----------------------------------------------------------------------------
x = x/Dims(3,1) + Dims(4,1);
y = y/Dims(3,2) + Dims(4,2);
z = z/Dims(3,3) + Dims(4,3);

X = x'*ones(1,V.dim(2));
Y = ones(V.dim(1),1)*y;

bbX = spm_dctmtx(Dims(1,1),Dims(2,1),x-1);
bbY = spm_dctmtx(Dims(1,2),Dims(2,2),y-1);
bbZ = spm_dctmtx(Dims(1,3),Dims(2,3),z-1);

dbX = spm_dctmtx(Dims(1,1),Dims(2,1),x-1,'diff');
dbY = spm_dctmtx(Dims(1,2),Dims(2,2),y-1,'diff');
dbZ = spm_dctmtx(Dims(1,3),Dims(2,3),z-1,'diff');

% Start progress plot
%----------------------------------------------------------------------------
spm_progress_bar('Init',length(z),'Modulating by determinants of Jacobian','planes completed');

M   = MF*Affine*inv(MG);
sgn = sign(det(M(1:3,1:3)));

img = zeros(V.dim(1:3));
% Cycle over planes
%----------------------------------------------------------------------------
for j=1:length(z)

	% Nonlinear deformations
	%----------------------------------------------------------------------------
	% 2D transforms for each plane
	tbx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ(j,:)', Dims(2,1), Dims(2,2) );
	tby = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ(j,:)', Dims(2,1), Dims(2,2) );
	tbz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *bbZ(j,:)', Dims(2,1), Dims(2,2) );

	tdx = reshape( reshape(Transform(:,1),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ(j,:)', Dims(2,1), Dims(2,2) );
	tdy = reshape( reshape(Transform(:,2),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ(j,:)', Dims(2,1), Dims(2,2) );
	tdz = reshape( reshape(Transform(:,3),Dims(2,1)*Dims(2,2),Dims(2,3)) *dbZ(j,:)', Dims(2,1), Dims(2,2) );

	% Jacobian of transformation from template
	% to affine registered image.
	%---------------------------------------------
	j11 = dbX*tbx*bbY' + 1;
	j12 = bbX*tbx*dbY';
	j13 = bbX*tdx*bbY';

	j21 = dbX*tby*bbY';
	j22 = bbX*tby*dbY' + 1;
	j23 = bbX*tdy*bbY';

	j31 = dbX*tbz*bbY';
	j32 = bbX*tbz*dbY';
	j33 = bbX*tdz*bbY' + 1;

	% Combine Jacobian of transformation from
	% template to affine registered image, with
	% Jacobian of transformation from affine
	% registered image to original image.
	%---------------------------------------------
	J11 = M(1,1)*j11 + M(1,2)*j21 + M(1,3)*j31;
	J12 = M(1,1)*j12 + M(1,2)*j22 + M(1,3)*j32;
	J13 = M(1,1)*j13 + M(1,2)*j23 + M(1,3)*j33;

	J21 = M(2,1)*j11 + M(2,2)*j21 + M(2,3)*j31;
	J22 = M(2,1)*j12 + M(2,2)*j22 + M(2,3)*j32;
	J23 = M(2,1)*j13 + M(2,2)*j23 + M(2,3)*j33;

	J31 = M(3,1)*j11 + M(3,2)*j21 + M(3,3)*j31;
	J32 = M(3,1)*j12 + M(3,2)*j22 + M(3,3)*j32;
	J33 = M(3,1)*j13 + M(3,2)*j23 + M(3,3)*j33;

	% The determinant of the Jacobian reflects relative volume changes.
	%------------------------------------------------------------------
	detJ = J11.*(J22.*J33 - J23.*J32) - J21.*(J12.*J33 - J13.*J32) + J31.*(J12.*J23 - J13.*J22);

	img(:,:,j) = spm_slice_vol(V,spm_matrix([0 0 j]),V.dim(1:2),1).*detJ*sgn;
	%img(:,:,j) = detJ*sgn;

	spm_progress_bar('Set',j);
end

% Write the data
%------------------------------------------------------------------
VO = V;
[pth,nm,xt,vr] = fileparts(deblank(V.fname));
VO.fname       = fullfile(pth,['m' nm xt vr]);
if isfield(VO,'descrip'),
	VO.descrip = [VO.descrip ' - Jacobian modulated'];
else,
	VO.descrip = ['spm - Jacobian modulated'];
end;
VO = spm_write_vol(VO, img);
spm_progress_bar('Clear');
return;
%_______________________________________________________________________