readSBItrials

 function [samp,pick1,pick2,bexp,ier]=readSBItrials(sqddat,sysid,trials,chgdcm)


 readSBI : Rev.1.0, 2001-11-17
 read SBI data
 function [samp,pick1,pick2,bexp,ier]=readSBI(sqddat,sysid,chgdcm)
 Input: sqddat  % $SQDDAT  ex. /data1/toyama2/
        sysid   % System ID  ex. 1000987
        chgdcm  % if 1, pick1/2 are converted to DICOM cordinate.
 Output: samp struct(nch,       % number of active gardiometers
                     nsamp,     % sampling number of measurement
                     nprsamp,   % sampling number for pre-trigger
                     sampf,     % sampling frequency
                     nrept,     % repeat number
                     sph0,      % origin of sphere model on DICOM cord.
                     sphr,      % radius of sphere model [m]
                     mrisys)    % MRI system ID
         pick1(ch,samp)         % position and direction of detection coils
                                   on DICOM cord.  [m],[normal unit vector]
         pick2(ch,samp)         % position and direction of compensation coils
                                   on DICOM cord.  [m],[normal unit vector]
         bexp(ch,samp,rept)  % measured fields  [T]
 history
 2001-11-09 S.Kajihara

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

0001 function [samp,pick12,pick22,bexp,ier]=readSBItrials(meg_dir,sysid,trials,chgdcm)
0002 % function [samp,pick1,pick2,bexp,ier]=readSBItrials(sqddat,sysid,trials,chgdcm)
0003 %
0004 %
0005 % readSBI : Rev.1.0, 2001-11-17
0006 % read SBI data
0007 % function [samp,pick1,pick2,bexp,ier]=readSBI(sqddat,sysid,chgdcm)
0008 % Input: sqddat  % $SQDDAT  ex. /data1/toyama2/
0009 %        sysid   % System ID  ex. 1000987
0010 %        chgdcm  % if 1, pick1/2 are converted to DICOM cordinate.
0011 % Output: samp struct(nch,       % number of active gardiometers
0012 %                     nsamp,     % sampling number of measurement
0013 %                     nprsamp,   % sampling number for pre-trigger
0014 %                     sampf,     % sampling frequency
0015 %                     nrept,     % repeat number
0016 %                     sph0,      % origin of sphere model on DICOM cord.
0017 %                     sphr,      % radius of sphere model [m]
0018 %                     mrisys)    % MRI system ID
0019 %         pick1(ch,samp)         % position and direction of detection coils
0020 %                                   on DICOM cord.  [m],[normal unit vector]
0021 %         pick2(ch,samp)         % position and direction of compensation coils
0022 %                                   on DICOM cord.  [m],[normal unit vector]
0023 %         bexp(ch,samp,rept)  % measured fields  [T]
0024 % history
0025 % 2001-11-09 S.Kajihara
0026 %
0027 % Copyright (C) 2011, ATR All Rights Reserved.
0028 % License : New BSD License(see VBMEG_LICENSE.txt)
0029 
0030 global vbmeg_inst;
0031 define = vbmeg_inst.const;
0032 
0033 pos_dir = '../pos';     % MEGディレクトリからposディレクトリへの相対path
0034 head_dir = '../head';   % MEGディレクトリからheadディレクトリへの相対path
0035 mri_dir = '../mri';
0036 
0037 ier=0;
0038 
0039 %
0040 % --- read SBI header
0041 %
0042 % file info.
0043 fid=fopen(sprintf('%s%s%s%s%s',meg_dir,head_dir,filesep,sysid,...
0044           define.HEADER_EXTENSION),'r','l');
0045 
0046 if fid==-1, disp('XXX Error : Cannot open head file !!\n'); ier=1; return; end
0047 
0048 H_file=meg_hread('file',fid);
0049 
0050 % condition
0051 H_meg=meg_hread('meg',fid,H_file.i_meg_begin);
0052 nch0=H_meg.i_nch;
0053 nsamp=H_meg.i_nsample;  % sampling number of measurement
0054 nprsamp=H_meg.i_nprsample;  % sampling number for pre-trigger
0055 sampf=H_meg.i_samp_freq;  % sampling frequency
0056 nrept=H_meg.i_nwave;  % repeat number
0057 vbit=H_meg.f_ad_weight;
0058 
0059 % all sensors
0060 kind0=[]; lock0=[]; pick10=[]; pick20=[]; amp0=[]; sens0=[];
0061 
0062 for i=1:nch0
0063     H_meg_ch=meg_hread('meg_ch',fid,H_file.i_meg_ch_begin+H_file.i_meg_ch_size*(i-1));
0064     kind0(i)=H_meg_ch.i_pick_type; lock0(i)=H_meg_ch.i_meg_lock;
0065     pick10(i,1:3)=(H_meg_ch.st_dcp)'; pick10(i,4:6)=(H_meg_ch.st_dca)';
0066     pick20(i,1:3)=(H_meg_ch.st_ccp)'; pick20(i,4:6)=(H_meg_ch.st_cca)';
0067     amp0(i)=H_meg_ch.f_amp_gain; sens0(i)=H_meg_ch.f_phy_v_coef;
0068 end
0069 
0070 if chgdcm,
0071     % positioning
0072     H_pos=meg_hread('pos',fid,H_file.i_pos_begin);
0073     npos=H_pos.i_ncoil;  % number of positioning coils
0074     H_pos_rept=meg_hread('pos_rept',fid,H_file.i_pos_rept_begin);
0075     pos_lcz_begin=H_pos_rept.i_data_begin;
0076     % MRI
0077     H_mri=meg_hread('mri',fid,H_file.i_mri_begin);
0078     mark=H_mri.st_axes(1:npos,:);  % MRI marker position on DICOM-shift cord.
0079     mrisys=H_mri.c_sysid{1};  % MRI system ID
0080     % LCZ
0081     H_lcz=meg_hread('lcz',fid,H_file.i_lcz_begin);
0082     sphr=H_lcz.f_model_r;  % radias of sphere model
0083     sph0=H_lcz.st_model;  % origin of sphere model on DICOM-shift cord.
0084 end
0085 
0086 fclose(fid); clear H_file H_meg H_meg_ch H_pos H_mri H_lcz;
0087 
0088 %
0089 % --- read SBI data
0090 %
0091 % meg data
0092 fid=fopen(sprintf('%s%s%s',meg_dir,sysid,...
0093           define.MEG1_EXTENSION),'r','l');
0094 
0095 if fid==-1, disp('XXX Error : Cannot open meg file !!\n'); ier=1; return; end
0096 
0097 a=fread(fid,nrept*nch0*nsamp,'short');
0098 fclose(fid);
0099 
0100 bexp00=(reshape(a,nch0*nsamp,nrept)); clear a;
0101 bexp0=bexp00(:,trials); clear bexp00;
0102 nrept=length(trials);
0103 bexp0=(reshape(bexp0,nch0,nsamp,nrept)); clear a;
0104 % -- patch to reject offset
0105 if 0,
0106     for i=1:nrept,
0107         for j=1:nch0,
0108             bexp0(j,:,i)=bexp0(j,:,i)-mean(bexp0(j,1:100,i));
0109         end
0110     end
0111 end
0112 
0113 if chgdcm,
0114     % pos. data
0115     fid=fopen(sprintf('%s%s%s%s%s',meg_dir,pos_dir,filesep,sysid,...
0116               define.POS_EXTENSION),'r','l');
0117 
0118     if fid==-1, 
0119         disp('XXX Error : Cannot open pos file !!\n'); ier=1; 
0120         return; 
0121     end
0122     
0123     fseek(fid,pos_lcz_begin,'bof');
0124     pos=[];
0125     for i=1:npos
0126         a=fread(fid,8,'float');
0127         b=reshape(a,1,8);
0128         pos(i,:)=b(1:3);  % position of positioning coils on DEWAR cord.
0129     end
0130     fclose(fid); clear a b;
0131 end
0132 
0133 %
0134 % ---  select active coils and data
0135 %
0136 okix=find(kind0==1 & lock0==0);
0137 nch=size(okix,2);  % number of active gardiometers
0138 pick1=pick10(okix,:);  % position and direction of detection coils on DEWAR cord. ([m],[rad,])
0139 pick2=pick20(okix,:);  % position and direction of compensation coils on DEWAR cord. ([m],[rad,])
0140 
0141 for i=1:nrept,
0142     bexp(:,:,i)=bexp0(okix,:,i)./(repmat(amp0(okix),nsamp,1))'.*(repmat(sens0(okix),nsamp,1))'*vbit;  % fields
0143 end
0144 
0145 clear kind0 lock0 pick10 pick20 amp0 sens0 vbit bexp0 okix;
0146 
0147 if ~chgdcm,
0148     pick12=pick1; pick22=pick2;
0149     samp=struct('nch',nch, ...  % number of active gardiometers
0150         'nsamp',nsamp, ...      % sampling number of measurement
0151         'nprsamp',nprsamp, ...  % sampling number for pre-trigger
0152         'sampf',sampf, ...      % sampling frequency
0153         'nrept',nrept);         % repeat number
0154     return;
0155 end
0156 
0157 %
0158 % --- convert cordinate from DICOM-shift or Dewar to DICOM
0159 %
0160 
0161 % origin of DICOM-shift
0162 fname1=sprintf('%s%s%s%s%s%s%s',meg_dir,mri_dir,filesep,mrisys,...
0163            filesep,mrisys,define.DICOM1_EXTENSION);
0164 
0165 fid=fopen(fname1);
0166 
0167 if fid==-1, 
0168     disp('XXX Error : Cannot open mri-table file !!\n'); ier=1; 
0169     return; 
0170 end
0171 
0172 a=fscanf(fid,'%s'); 
0173 fclose(fid);
0174 
0175 sx=str2num(a(findstr(a,'image_x_vector=')+15:findstr(a,'image_y_vector=')-1));
0176 sy=str2num(a(findstr(a,'image_y_vector=')+15:findstr(a,'image_z_vector=')-1));
0177 sz=str2num(a(findstr(a,'image_z_vector=')+15:findstr(a,'pixel_spacing=')-1));
0178 s0=str2num(a(findstr(a,'z0_slice_origin=')+16:findstr(a,'magnification=')-1));
0179 ss=str2num(a(findstr(a,'pixel_spacing=')+14:findstr(a,'z0_slice_origin=')-1));
0180 sp1=findstr(a,'slice_pitch=')+12; sp2=findstr(a(sp1:size(a,2)),'te=')-1; sp=str2num(a(sp1:sp1+sp2-1));
0181 nsl=str2num(a(findstr(a,'number_of_slices=')+17:findstr(a,'number_of_echos=')-1));
0182 
0183 b=1; 
0184 clear sloc; 
0185 fid=fopen(fname1);
0186 
0187 while b,
0188     if ~isempty(findstr(char(fgets(fid)),'slice_location_table')), b=0; end
0189 end
0190 
0191 for i=1:nsl, sloc(i)=str2num(char(fgets(fid))); end
0192 fclose(fid);
0193 scorner=sz*min(sloc)+s0;
0194 dcms0=128*sx*ss+127*sy*ss+fix(nsl/2)*sz*sp+scorner; dcms0=dcms0.*1e-3;
0195 % from DICOM-shift to DICOM
0196 mark2=mark+repmat(dcms0,npos,1);  % MRI marker position on DICOM cord.
0197 if sphr<0, sph02=sph0+dcms0';  % origin of sphere model on DICOM cord.
0198 else, sph02=[0,0,0]; disp('XXX : No suport, because not using VIVID'); ier=1; return; end
0199 % from Dewar to DICOM
0200 pick1(:,4)=mod(pick1(:,4),pi); pick1(:,5:6)=mod(pick1(:,5:6),2*pi); % normalize 0<pi or 2pi
0201 pick2(:,4)=mod(pick2(:,4),pi); pick2(:,5:6)=mod(pick2(:,5:6),2*pi); % normalize 0<pi or 2pi
0202 pick1(:,4:6)=[sin(pick1(:,4)).*cos(pick1(:,5)) sin(pick1(:,4)).*sin(pick1(:,5)) cos(pick1(:,4))]; % angle -> normal vector
0203 pick2(:,4:6)=[sin(pick2(:,4)).*cos(pick2(:,5)) sin(pick2(:,4)).*sin(pick2(:,5)) cos(pick2(:,4))]; % angle -> normal vector
0204 for i=1:npos  % minimum distance between a positioning coil and detection coils
0205     sfromc(i)=min(sqrt(sum((pick1(:,1:3)-repmat(pos(i,:),size(pick1,1),1)).^2,2)));
0206 end
0207 sel=find(sfromc<max(sfromc)); % select 3 coils
0208 ex=pos(sel(3),:)-pos(sel(1),:); ex=ex/norm(ex);
0209 ez=cross(ex,pos(sel(2),:)-pos(sel(1),:)); ez=ez/norm(ez);
0210 ey=cross(ez,ex); mc=[ex; ey; ez]; % axes of Head cord. on Dewar cord.
0211 oc=sum((pos(sel(2),:)-pos(sel(1),:)).*ex)*ex + pos(sel(1),:); % origin of Head cord. on Dewar cord.
0212 ex=mark2(sel(3),:)-mark2(sel(1),:); ex=ex/norm(ex);
0213 ez=cross(ex,mark2(sel(2),:)-mark2(sel(1),:)); ez=ez/norm(ez);
0214 ey=cross(ez,ex); mm=[ex; ey; ez]; % axes of Head cord. on DICOM cord.
0215 om=sum((mark2(sel(2),:)-mark2(sel(1),:)).*ex)*ex + mark2(sel(1),:); % origin of Head cord. on DICOM cord.
0216 trans = mm'*mc;
0217 pick12(:,1:3)=(trans*(pick1(:,1:3)-repmat(oc,nch,1))')'+repmat(om,nch,1);  % position of detection coils on DICOM cord.
0218 pick22(:,1:3)=(trans*(pick2(:,1:3)-repmat(oc,nch,1))')'+repmat(om,nch,1);  % position of compensation coils on DICOM cord.
0219 pick12(:,4:6)=(trans*pick1(:,4:6)')'; % direction of detection coils on DICOM cord.
0220 pick22(:,4:6)=(trans*pick2(:,4:6)')'; % direction of compensation coils on DICOM cord.
0221 samp=struct('nch',nch, ...  % number of active gardiometers
0222     'nsamp',nsamp, ...      % sampling number of measurement
0223     'nprsamp',nprsamp, ...  % sampling number for pre-trigger
0224     'sampf',sampf, ...      % sampling frequency
0225     'nrept',nrept, ...      % repeat number
0226     'sph0',sph02, ...       % origin of sphere model on DICOM cord.
0227     'sphr',sphr, ...        % radius of sphere model [m]
0228     'mrisys',mrisys);       % MRI system ID
0229