readSBI_yosio

 readSBI : Rev.1.0, 2001-11-17
 read SBI data
 function [samp,pick1,pick2,bexp,ier]=readSBI(meg_dir,sysid,chgdcm)
 Input: meg_dir  % $MEG_DIR  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

 Modified by Taku Yoshioka
 2003-10-09 
 ・ファイル拡張子を構造体defineから取得

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

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