function M = spm_get_space5(P,M)
% Get/set the voxel-to-world mapping of an image
% FORMAT M = spm_get_space(P)
%            spm_get_space(P,M)
% M - voxel-to-world mapping
% P - image filename
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience

% John Ashburner
% $Id: spm_get_space5.m,v 1.1.1.1 2006/01/09 19:29:21 kpearson Exp $


[pth,nam,ext] = fileparts(P);
t = find(ext==',');
n = [1 1];
if ~isempty(t),
    n   = str2num(ext((t(1)+1):end));
    ext = ext(1:(t-1));
    P   = fullfile(pth,[nam ext]);
end;

N = nifti(P);


if (isempty(t) | n == 1)  & (strcmpi(ext,'.img') | strcmpi(ext,'.img') )
	aux_file = deblank(N.aux_file);
	if strcmp(aux_file,'none'),
		v = spm_vol_nifti(P); 
		[normalized,wfu_mat,wfu_M,spm99_norm] = wfu_isnormalized(P,v);
		disp('spm_get_space SPM99 file detected');
		N.mat = wfu_mat;
		N.mat0 = wfu_mat;
	end
end
if nargin==2,
    N.mat_intent = 'Aligned';
    if n(1)==1,

        % Ensure volumes 2..N have the original matrix
        if size(N.dat,4)>1 && sum(sum((N.mat-M).^2))>1e-8,
            M0 = N.mat;
            if ~isfield(N.extras,'mat'),
                N.extras.mat = zeros([4 4 size(N.dat,4)]);
            else
                if size(N.extras.mat,4)<size(N.dat,4),
                    N.extras.mat(:,:,size(N.dat,4)) = zeros(4);
                end;
            end;
            for i=2:size(N.dat,4),
                if sum(sum(N.extras.mat(:,:,i).^2))==0,
                    N.extras.mat(:,:,i) = M0;
                end;
            end;
        end;

        N.mat        = M;
        if strcmp(N.mat0_intent,'Aligned'), N.mat0 = M; end;
        if ~isempty(N.extras) && isstruct(N.extras) && isfield(N.extras,'mat') &&...
            size(N.extras.mat,3)>=1,
            N.extras.mat(:,:,n(1)) = M;
        end;
    else
        N.extras.mat(:,:,n(1)) = M;
    end;
    create(N);
else
    if ~isempty(N.extras) && isstruct(N.extras) && isfield(N.extras,'mat') &&...
        size(N.extras.mat,3)>=n(1) && sum(sum(N.extras.mat(:,:,n(1)).^2)),
        M = N.extras.mat(:,:,n(1));
    else
        M = N.mat;
    end;
end;

