readSBI

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