vb_BVMake2D3DTRF

 function vb_BVMake2D3DTRF(TRFfn,T1fn,EPIfn,VMRdim)

 coregister & make 2D-3D TRF file
 (this function dos not modify original mat files)

 TRFfn: TRF filename (*.trf)
  T1fn: T1 3D filename (*.img)
  EPIfn: EPI filename (*.img)
  VMRdim(option): VMR dimension ([R-L, A-P, S-I]), default [256 256 256]


 2002/12/19 N.Goda

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

0001 function vb_BVMake2D3DTRF(TRFfn,T1fn,EPIfn,VMRdim)
0002 % function vb_BVMake2D3DTRF(TRFfn,T1fn,EPIfn,VMRdim)
0003 %
0004 % coregister & make 2D-3D TRF file
0005 % (this function dos not modify original mat files)
0006 %
0007 % TRFfn: TRF filename (*.trf)
0008 %  T1fn: T1 3D filename (*.img)
0009 %  EPIfn: EPI filename (*.img)
0010 %  VMRdim(option): VMR dimension ([R-L, A-P, S-I]), default [256 256 256]
0011 %
0012 %
0013 % 2002/12/19 N.Goda
0014 %
0015 % Copyright (C) 2011, ATR All Rights Reserved.
0016 % License : New BSD License(see VBMEG_LICENSE.txt)
0017 
0018 if nargin<3
0019   fprintf('\n invalid arguments\n');
0020   return
0021 elseif nargin<4
0022   VMRdim=[256 256 256];
0023 end
0024 
0025 % preparation -------------------------------------------------------------
0026 SWD=spm('Dir');
0027 global sptl_Ornt;
0028 if ~isempty(sptl_Ornt)
0029   tmpglobal=sptl_Ornt; % copy
0030   restoreglobal=1;
0031 else
0032   restoreglobal=0;
0033 end
0034 sptl_Ornt = [0 0 0  0 0 0 -1 1 1 0 0 0]; % Radiological Convention (L is R)
0035 DIR1 = [fullfile(SWD,'templates'),filesep];
0036 
0037 % copy files
0038 delete tmpt1* tmpepi*;
0039 % Standard T1 3D
0040 copyfile(T1fn,'tmpt1.img');
0041 copyfile([T1fn(1:end-3) 'hdr'],'tmpt1.hdr');
0042 T1F='tmpt1.img';
0043 T1T=fullfile(DIR1,'T1.img');
0044 % EPI
0045 copyfile(EPIfn,'tmpepi.img');
0046 copyfile([EPIfn(1:end-3) 'hdr'],'tmpepi.hdr');
0047 EPIF='tmpepi.img';
0048 EPIT=fullfile(DIR1,'EPI.img');
0049 
0050 % read header info
0051 % target (T1)
0052 [T1dim T1vox]=spm_hread(T1F);
0053 if ~isequal(T1vox,[1 1 1])
0054   fprintf('\n T1 3D should be 1mm^3 voxel\n');
0055   return
0056 end
0057 T1origin0=(T1dim(1:3)+1)/2; 
0058 T1originmm0=T1vox.*T1origin0; % ex. [96 128.5 128.5] in mm
0059 % object (EPI)
0060 [EPIdim EPIvox]=spm_hread(EPIF);
0061 EPIorigin0=(EPIdim(1:3)+1)/2; % dim64->32.5, % dim50->25.5
0062 EPIoriginmm0=EPIvox.*EPIorigin0;
0063 
0064 % default BV center (maybe) [0-255], ex. [95,128,128] if VMRdim==[191,256,256]
0065 BVcenter=round((VMRdim-1)/2);
0066 
0067 % default BV center in SPM-voxel Coordinate in target space [1-256]
0068 TGTV=BVvox2spmv(BVcenter,T1dim,VMRdim);   % ex. [96,128,128] if T1dim==[191,256,256]
0069 originmmoff=TGTV-T1originmm0; % origin offset in mm
0070 
0071 % set defalut mat file (adjust origin to BV origin) -------------------------
0072 T1originmm1=T1originmm0+originmmoff;
0073 off=-T1originmm1; % mm, origin adjusted
0074 MT1=[T1vox(1) 0 0 off(1) ; 0 T1vox(2) 0 off(2) ; 0 0 T1vox(3) off(3) ; 0 0 0 1];
0075 spm_get_space(T1F,MT1); % set default mat
0076 
0077 % EPI
0078 EPIoriginmm1=EPIoriginmm0+originmmoff;
0079 off=-EPIoriginmm1; % mm, origin adjusted
0080 M0=[EPIvox(1) 0 0 off(1) ; 0 EPIvox(2) 0 off(2) ; 0 0 EPIvox(3) off(3) ; 0 0 0 1];
0081 spm_get_space(EPIF,M0); % set default mat
0082 
0083 % coregister EPI -> T1 -----------------------------------------------------
0084 % pre-coregitered OBJ center in SPM-voxel Coordinate in object space
0085 OBJV0=inv(M0)*[0 0 0 1]'; 
0086 
0087 % first step: default coregister (rough coregistration)
0088 fprintf('\n first step coregistration ...');
0089 spm_coregister(T1F, EPIF, T1T, EPIT, [],'n');
0090 M1=spm_get_space(EPIF); % get coregistered mat
0091 
0092 % second step: MI coregister
0093 fprintf('\n MI coregistration ...');
0094 x=spm_mireg(spm_vol(T1F), spm_vol(EPIF));
0095 M2=inv(spm_matrix(x));
0096 
0097 % final result
0098 M3=M2*M1;
0099 spm_get_space(EPIF, M3); %  overwrite by new mat
0100 
0101 % convert to TRF -------------------------------------------------------------
0102 % coregistered OBJ center in SPM-mm Coordinate
0103 x1=M3*OBJV0;
0104 
0105 % convert to SPM-voxel Coordinates in target space (1-256)
0106 x2=inv(MT1)*x1;
0107 
0108 % convert to BV-VMR Coordinates (0-255)
0109 trf(1:3)=BVspmv2vox(x2(1:3)',T1dim,VMRdim);
0110 
0111 % rotation
0112 x=spm_imatrix(M3);
0113 trf(4)=x(4)/pi*180;
0114 trf(5)=-x(5)/pi*180;
0115 trf(6)=-x(6)/pi*180;
0116 
0117 % result output
0118 fid(1)=1;
0119 fid(2)=fopen(TRFfn,'w');
0120 fprintf('\nTRF\n');
0121 for i=1:2
0122   fprintf(fid(i),'\n');
0123   fprintf(fid(i),'y: %f\n',trf(2));
0124   fprintf(fid(i),'z: %f\n',trf(3));
0125   fprintf(fid(i),'x: %f\n',trf(1));
0126   fprintf(fid(i),'TraCor: %f\n',trf(4));
0127   fprintf(fid(i),'TraSag: %f\n',trf(5));
0128   fprintf(fid(i),'SagCor: %f\n',trf(6));
0129   fprintf(fid(i),'xFoV: %f\n',EPIdim(1)*EPIvox(1));
0130   fprintf(fid(i),'yFoV: %f\n',EPIdim(2)*EPIvox(2));
0131   fprintf(fid(i),'zFoV: %f\n',EPIdim(3)*EPIvox(3));
0132   fprintf(fid(i),'n_slices: %d\n',EPIdim(3));
0133   fprintf(fid(i),'sl_thick: %f\n',EPIvox(3));
0134   fprintf(fid(i),'sl_gap:   0.000000\n');
0135 end
0136 fclose(fid(2));
0137 
0138 if restoreglobal
0139   sptl_Ornt=tmpglobal; 
0140 end
0141 
0142 return