vb_set_head_model

 
 Set head and dipole model for MEG forward calculation test

 rdir = [getenv('MATHOME') '/SBIdata-new/neuroimage_retino/data/'];
 udir= [getenv('MATHOME') '/SBIdata-new/BEM/head/'];

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

0001 %
0002 % Set head and dipole model for MEG forward calculation test
0003 %
0004 % rdir = [getenv('MATHOME') '/SBIdata-new/neuroimage_retino/data/'];
0005 % udir= [getenv('MATHOME') '/SBIdata-new/BEM/head/'];
0006 %
0007 % Copyright (C) 2011, ATR All Rights Reserved.
0008 % License : New BSD License(see VBMEG_LICENSE.txt)
0009 
0010 rdir = 'C:\SBIdata-new\neuroimage_retino\data\';
0011 
0012 headfile  = '100008d.s01.mat';            
0013 brainfile = 'brain-retino2';
0014 pickfile  = 'meg-pick';
0015 
0016 if exist('Nsphere','var') & ~isempty(Nsphere)
0017     [Fhead,XXhead] = vb_make_fullsphere(Nsphere);
0018     Vhead = XXhead;
0019 else
0020     %load([udir headfile]);
0021     %[Fhead, Vhead, XXhead] = vb_out_normal(F,V);
0022 end
0023 
0024 load([rdir brainfile], 'V', 'xx');
0025 load([rdir pickfile], 'pick');
0026 
0027 if ~isempty(Ndipole)
0028     Nall = size(V,1);
0029     step = fix(Nall/Ndipole);
0030     V    = V(1:step:Nall,:);
0031     xx   = xx(1:step:Nall,:);
0032 end
0033 
0034 pick1 = [ pick(:,1:3)];    % センサー位置1
0035 pick2 = [ pick(:,4:6)];    % センサー位置2
0036 Qpick = [ pick(:,7:9)];    % センサー方向
0037 
0038 Npick    = size(pick,1);
0039 Ndipole = size(V,1);
0040 
0041 rpick = sum(sqrt(sum(pick1.^2,2)))/Npick;
0042 dpick = sum(sqrt(sum((pick1 - pick2).^2,2)))/Npick;
0043 
0044 rhead = 1.1;
0045 r     = sqrt(sum(V.^2,2));
0046 rmax  = max(r);
0047 R0    = rmax*rhead;
0048 
0049 if exist('Nsphere','var') & ~isempty(Nsphere)
0050     Vhead = R0*Vhead;
0051 end
0052 
0053 if EEGmode==0; return; end;
0054 
0055 %
0056 % 3層モデル
0057 %
0058 
0059 sigma = [ 1.0  0.0125  1.0  0.0];
0060 R     = [ 0.87  0.92   1.0 ];
0061 R     = R/R(1);
0062 NV    = size(Vhead,1);
0063 NF    = size(Fhead,1);
0064 
0065 Vhead0    = Vhead;
0066 Fhead0    = Fhead;
0067 XXhead0 = XXhead;
0068 Nvertex = [1 NV];
0069 Npatch    = [ 1 NF];
0070 
0071 for n=2:Nsurf
0072     Vhead    = [Vhead ; Vhead0 * R(n)];
0073     Fhead    = [Fhead ; Fhead0 + NV*(n-1)];
0074     XXhead    = [XXhead; XXhead0];
0075     Nvertex = [Nvertex; NV*(n-1)+1, NV*n];
0076     Npatch    = [Npatch;  NF*(n-1)+1, NF*n];
0077 end
0078 
0079 R            = R(1:Nsurf);
0080 BEM.R        = R;
0081 BEM.sigma    = sigma;
0082 BEM.Nvertex = Nvertex;
0083 BEM.Npatch  = Npatch;
0084 
0085 fprintf('# of surface  = %d\n',Nsurf);
0086 fprintf('# of vertices = %d\n',size(Vhead,1));
0087 
0088 if Nsurf==1; return; end;
0089 
0090 % EEG sensor
0091 IDsensor    = IDsensor + NV*(Nsurf-1);
0092 Xeeg        = Vhead(IDsensor,:);