 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;
 %_______________________________________________________________________