
% FORMAT   Yo = normalize_voxels
% __________________________________________________________________________
% This function takes in a set of *.img images, and produces the same number 
% of output image-structures, by dividing each voxel-value by _it's_ mean
% intensity[i.e., computed for that voxel over all input images]. These image 
% intensities are stored in a 4-D matrix [the 4th dimension indexes the images]
% called 'Yo' in the matlab workspace, as well as an spm-Vol structure, 'Vo'.
% [will soon modify code to save them as *.img files...]
% The input images must have the same dimensions, orientations (as defined 
% by the origin header field or any associated *.mat files), and the same
% voxel sizes.
%
% __________________________________________________________________________
%   Pavithra Srinivasan, 12th Sep 2002.
%

clear all;

%-Select images & check dimensions, orientations and voxel sizes
%-----------------------------------------------------------------------
fprintf('\t...select files')
P = spm_get(Inf,'.img','Select images to be normalized');
fprintf(' ...mapping & checking files')

Vi    = spm_vol(P);
nvox  = prod(size(Vi));%-total number of voxels in any one image.
nimgs = size(P, 1);    %-total number of selected images.

%-error checking:
if nvox==0, fprintf('\t%s : no images selected\n\n',mfilename), return, end
if nvox>1 & any(any(diff(cat(1,Vi.dim),1,1),1)&[1,1,1,0])
	error('images don''t all have same dimensions'), end
if any(any(any(diff(cat(3,Vi.mat),1,3),3)))
	error('images don''t all have same orientation & voxel size'), end

%-Reference....
%-----------------------------------------------------------------------
%Vmean = struct(	'fname',	'mean.img',...
%		'dim',		[Vi(1).dim(1:3),4],...
%		'mat',		Vi(1).mat,...
%		'pinfo',	[1.0,0,0]',...
%		'descrip',	'spm - mean image');


[Yi]     = spm_read_vols(Vi);
Ymean    = mean(Yi, 4); %- Ymean is now the mean 3-D intensity-matrix.
Vo       = Vi;

%-Dividing each voxel by "it's" mean:
for i = 1 : nimgs
    Yo(:, :, :, i) = Yi(:, :, :, i) ./ Ymean;
end

% The result 'Yo' which is a 4-D matrix[the 4th dimension
% indexes the images] contains the new normalized intensity values.

%-Creating image files from the Yo matrix:
resdir = spm_get(-1,'.','Select directory to save output files.')
for i = 1 : nimgs
    Vo(i).fname = prepend(Vo(i).fname, 'voxnorm_', resdir);
    Vo(i)       = spm_write_vol(Vo(i), Yo(:, :, :, i));
end

%clear Vi; 
%clear Yi;
fprintf(' ...done\n')



