vb_unsn

 ---
 function [XYZmm]=vb_unsn(nXYZmm,sn)

 nXYZmm: mm coordinates in normalized volume [ 3 x nvoxels ]
 sn: filename or structure

 XYZmm: mm coordinates in original volume    [ 3 x nvoxels ]

 2004/12/14 N.Goda
 2006/02/23 Kawawaki, Dai
 2008/09/11 Kawawaki, Dai
 2010/05/12 Y. Takeda added spm8 support
 2010/07/22 M. Sato simplified the switch logic, XYZvx is removed

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function [XYZmm]=vb_unsn(nXYZmm,sn)
0002 % ---
0003 % function [XYZmm]=vb_unsn(nXYZmm,sn)
0004 %
0005 % nXYZmm: mm coordinates in normalized volume [ 3 x nvoxels ]
0006 % sn: filename or structure
0007 %
0008 % XYZmm: mm coordinates in original volume    [ 3 x nvoxels ]
0009 %
0010 % 2004/12/14 N.Goda
0011 % 2006/02/23 Kawawaki, Dai
0012 % 2008/09/11 Kawawaki, Dai
0013 % 2010/05/12 Y. Takeda added spm8 support
0014 % 2010/07/22 M. Sato simplified the switch logic, XYZvx is removed
0015 %
0016 % Copyright (C) 2011, ATR All Rights Reserved.
0017 % License : New BSD License(see VBMEG_LICENSE.txt)
0018 
0019 if isstr(sn), % backward compatibility
0020   sntmp=load(sn);
0021   sn=sntmp; 
0022 end
0023 if isfield(sn,'Dims')
0024   ver = 'spm99';
0025 else
0026   switch lower(sn.VG(1).descrip)
0027    case 'minc file'
0028     ver = 'spm2';
0029    case 'nifti-1 image'
0030     ver = 'spm5';
0031    case 'grey matter tissue probability'
0032     ver = 'spm5';
0033    otherwise
0034     error('Cannot get SPM version.');
0035   end
0036 end
0037 
0038 switch ver
0039  case 'spm99'
0040   tmpl_dim=sn.Dims(1,:);    % [91,109,91]
0041   tmpl_nbasis=sn.Dims(2,:);
0042   tmpl_vox=sn.Dims(3,:);    % [2 2 2]
0043   tmpl_origin=sn.Dims(4,:); % [46,64,37]
0044   
0045   org_dim= sn.Dims(5,:);
0046   org_vox=sn.Dims(6,:);
0047 
0048   MF        = sn.MF;
0049   Transform = reshape(sn.Transform,tmpl_nbasis(1),...
0050               tmpl_nbasis(2),tmpl_nbasis(3),3); % 2d->4d
0051   Affine    = sn.Affine;
0052 
0053  case {'spm2'}
0054   tmpl_dim=sn.VG(1).dim(1:3); % [91,109,91]
0055   tmpl_nbasis=size(sn.Tr);
0056   tmpl_vox=sqrt(sum(sn.VG(1).mat(1:3,1:3).^2));  % [2 2 2]
0057   tmpl_origin=sn.VG(1).mat\[0 0 0 1]'; %[46,64,37]
0058   tmpl_origin=tmpl_origin(1:3)';
0059 
0060   org_dim=sn.VF.dim(1:3);  % [191,256,256] etc
0061   org_vox=sqrt(sum(sn.VF.mat(1:3,1:3).^2));  % [2 2 2]
0062 
0063   MF        = sn.VF.mat;
0064   Transform = sn.Tr; % 4d
0065   Affine    = sn.Affine;
0066  case {'spm5'}
0067      XYZmm = spm_get_orig_coord_spm5(nXYZmm', sn)';
0068      return;
0069  otherwise
0070   error('There was a significant error. This program was aborted.');
0071 end
0072 
0073 flip_flag = 0;
0074 
0075 switch ver
0076  case {'spm5'}
0077   [U,S,V] = svd(Affine(1:3,1:3));
0078   if det(U*V') > 0
0079     flip_flag = 1;
0080   end
0081 end
0082 
0083 % voxel coordinate in template
0084 % (2mm^3 resolution, 91x109x91, origin [46,64,37])
0085 % non-integer values are allowed
0086 x=nXYZmm(1,:)/tmpl_vox(1)+tmpl_origin(1);
0087 y=nXYZmm(2,:)/tmpl_vox(2)+tmpl_origin(2);
0088 z=nXYZmm(3,:)/tmpl_vox(3)+tmpl_origin(3);
0089 
0090 basX = spm_dctmtx(tmpl_dim(1),tmpl_nbasis(1),x-1);
0091 basY = spm_dctmtx(tmpl_dim(2),tmpl_nbasis(2),y-1);
0092 basZ = spm_dctmtx(tmpl_dim(3),tmpl_nbasis(3),z-1);
0093 
0094 Mult = MF*Affine;
0095 XYZmm=zeros(size(nXYZmm));
0096 
0097 trx0=reshape(Transform(:,:,:,1),tmpl_nbasis(1)*tmpl_nbasis(2),...
0098          tmpl_nbasis(3));
0099 try0=reshape(Transform(:,:,:,2),tmpl_nbasis(1)*tmpl_nbasis(2),...
0100          tmpl_nbasis(3));
0101 trz0=reshape(Transform(:,:,:,3),tmpl_nbasis(1)*tmpl_nbasis(2),...
0102          tmpl_nbasis(3));
0103 
0104 for j=1:length(z)
0105   tx0 = trx0*basZ(j,:)';
0106   ty0 = try0*basZ(j,:)';
0107   tz0 = trz0*basZ(j,:)';
0108   tx1 = reshape(tx0,tmpl_nbasis(1),tmpl_nbasis(2)) *basY(j,:)';
0109   ty1 = reshape(ty0,tmpl_nbasis(1),tmpl_nbasis(2)) *basY(j,:)';
0110   tz1 = reshape(tz0,tmpl_nbasis(1),tmpl_nbasis(2)) *basY(j,:)';
0111   tx2 = tx1' *basX(j,:)';
0112   ty2 = ty1' *basX(j,:)';
0113   tz2 = tz1' *basX(j,:)';
0114 
0115   X1 = x(j)+tx2;
0116   Y1 = y(j)+ty2;
0117   Z1 = z(j)+tz2;
0118 
0119   X2= Mult(1,1)*X1 + Mult(1,2)*Y1 + Mult(1,3)*Z1 + Mult(1,4);
0120   Y2= Mult(2,1)*X1 + Mult(2,2)*Y1 + Mult(2,3)*Z1 + Mult(2,4);
0121   Z2= Mult(3,1)*X1 + Mult(3,2)*Y1 + Mult(3,3)*Z1 + Mult(3,4);
0122   
0123   XYZmm(:,j)=[X2,Y2,Z2]';
0124   
0125   if flip_flag,
0126     XYZmm(1,j) = -XYZmm(1,j);
0127   end
0128 
0129 end
0130 
0131 return